A slow-down time-transformed symplectic integrator for solving the few-body problem

Long Wang^1,2, Keigo Nitadori² and Junichiro Makino²
¹Department of Astronomy, School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo, 113-0033, Japan
²RIKEN Center for Computational Science, 7-1-26 Minatojima-minami-machi, Chuo-ku, Kobe, Hyogo 650-0047, Japan
E-mail:[email protected]

(Accepted – . Received –; in original form –)

Abstract

An accurate and efficient method dealing with the few-body dynamics is important for simulating collisional $N$ -body systems like star clusters and to follow the formation and evolution of compact binaries. We describe such a method which combines the time-transformed explicit symplectic integrator (Preto & Tremaine, 1999; Mikkola & Tanikawa, 1999) and the slow-down method (Mikkola & Aarseth, 1996). The former conserves the Hamiltonian and the angular momentum for a long-term evolution, while the latter significantly reduces the computational cost for a weakly perturbed binary. In this work, the Hamilton equations of this algorithm are analyzed in detail. We mathematically and numerically show that it can correctly reproduce the secular evolution like the orbit averaged method and also well conserve the angular momentum. For a weakly perturbed binary, the method is possible to provide a few order of magnitude faster performance than the classical algorithm. A publicly available code written in the c++ language, sdar, is available on GitHub. It can be used either as a standalone tool or a library to be plugged in other $N$ -body codes. The high precision of the floating point to $62$ digits is also supported.

keywords:

methods: numerical – software: simulations – Galaxy: globular clusters: general

^†^†pagerange: A slow-down time-transformed symplectic integrator for solving the few-body problem–A^†^†pubyear: 2002

1 Introduction

Few-body dynamics embedded in $N$ -body systems like star clusters is important for many research topics: the formation of high-velocity stars via the ejections by interacting with binaries or massive black holes (e.g. Poveda, Ruiz & Allen, 1967; Leonard & Duncan, 1988, 1990; Yu & Tremaine, 2003; Gualandris, Portegies Zwart & Sipior, 2005; Gvaramadze, Gualandris & Portegies Zwart, 2009; Fujii & Portegies Zwart, 2011); the mergers of binaries that produce exotic objects like blue stragglers (e.g. Bailyn, 1995; Davies, Piotto & de Angeli, 2004; Hurley, et al., 2005; Heggie & Giersz, 2008; Leigh, Sills & Knigge, 2011; Hypki & Giersz, 2013) and gravitational wave sources (e.g. Portegies Zwart & McMillan, 2000; O’Leary et al., 2006; Banerjee, Baumgardt & Kroupa, 2010; Antonini, et al., 2016; Askar, Szkudlarek, Gondek-Rosińska, Giersz & Bulik, 2017; Di Carlo, et al., 2019); and the energy generation source that controls the dynamical evolution of star clusters (e.g. Hénon, 1961; Heggie, 1975; Hills, 1975; Spitzer, 1987; Binney & Tremaine, 1987; Breen & Heggie, 2013; Wang, 2020).

However, few-body systems are also challenging to handle in the numerical simulations. Firstly, they can have a much shorter dynamical timescale than that of the host environment. For example, an open cluster with a few thousand stars has a typical dynamical (crossing) time of few Myrs, but a close binary can have a period of few days. It is very time consuming to accurately follow the lifetime of the cluster with a time resolution less than the binary period. Secondly, the numerical error introduced by the integrator can accumulate after many binary orbits and cause an artificial drift of the orbital parameters.

To avoid these issues, the $N$ -body codes for simulating star clusters, such as nbody6 (Aarseth, 2003), artificially “freeze” the orbits of weakly perturbed binaries and hierarchical systems, i.e., if the perturbation is below a threshold and the hierarchical systems satisfy a stability criterion, the internal motion of the systems are not evolved until the perturbation becomes strong. Such trick can significantly reduce the computational cost, but it ignores the long-term effect of weak perturbation. In particular, the evolution of angular momentum, and thus the eccentricity, can be completely wrong. Besides, whether the instability of hierarchical systems is properly treated depends on the quality of the stability criterion.

There are two approximate solutions that can properly handle the few-body systems and also avoid the time consuming calculation. The first method is using the orbit averaged Hamiltonian. In the case of a hierarchical triple, instead of tracing the positions and velocities of individual components by integrating the classical Newtonian equation of motion, the orbits of inner and outer binaries can be evolved by using the Hamilton equation written in the Delaunay’s elements. The short-period terms in the Hamiltonian can be eliminated (orbit averaged) by the Von Zeipel transformation (e.g. Naoz, et al., 2013; Naoz, 2016). For example, the orbit averaged method with a quadruple-level perturbation term is provided in Naoz, et al. (2013). Because the timescale of secular evolution is longer than the periods of binaries, such method can be very efficient to integrate the orbital evolution of a stable hierarchical triple.

Another method is called “slow-down” introduced by Mikkola & Aarseth (1996). For a perturbed binary, by artificially slowing down the orbital motion (scaling the time) while keeping the orbital parameters unchanged, the external perturbation is effectively enlarged. In such case, the effect of perturbation on one binary orbit can represent the average effect of several orbits. When the perturbation is weak, this method can properly approximate the secular evolution. The benefit is rich: it is easy to be implemented in a $N$ -body code; can be applied for a general few-body system with an arbitrary number of binaries and singles; is not only for a kepler binary, but also for any type of periodic orbits. In Mikkola & Aarseth (1996), a constant slow-down factor is tested for a binary with a post-Newtonian type of perturbation. However, it is not well discussed and tested how we can change the slow-down factor, even though it is crucial for the actual use of the slow-down method. Especially, this is important for dealing with the common few-body systems in a star cluster, such as triples with high-eccentric or hyperbolic orbits of perturbers. When the perturbation to a binary becomes smaller, we start from no slow-down to a larger slow-down factor. But how to change it properly is not well understood. Mikkola & Aarseth (1996) discuss a factor-of-two change, but such a way would cause the non-conservative change in the total energy of the system.

These two methods not only can reduce the computing cost, but also avoid the large accumulate error due to much less integration steps. But as the price of approximation, both methods lose the phase information of the orbits, thus the mean motion resonance cannot be correctly reproduced. They also have different disadvantages. The orbit averaged method was used to study a limited type of stable hierarchical systems, e.g. stable hierarchical triples and the Lunar system. For a more general case like the hyperbolic encounters, systems with more than three bodies, or strongly perturbed few-body systems, to derive the equation of motion with high-order perturbation terms is unpractical. The slow-down method is much more flexible but it cannot avoid the issue of the singularity of Newtonian force, i.e., to integrate the motion of a high-eccentric binary, a large numerical error appears at the peri-center unless the integration step is significantly reduced.

The solution for the issue of singularity is to apply the regularization algorithm. The key idea of regularization is to remove the singularity by a transformation of the equation of motion. The Burdet-Heggie regularization (Burdet, 1967, 1968; Heggie, 1973) for solving the Kepler problem uses a time-transformation, $\mathrm{d}t=r\mathrm{d}\tau$ , where $\tau$ is a fictitious time. By using the eccentric vector in the equation of motion together, the singularity term can be eliminated. Kustaanheimo & Stiefel (1965) introduce the KS regularization method, which transforms the equation of motion of the Kepler orbit to a form of harmonic oscillator without singularity. Mikkola & Aarseth (1996) show how to combine the KS regularization together with the slow-down method to have an efficient solution for integrating a few-body system. Mikkola & Tanikawa (1999) and Preto & Tremaine (1999) introduce the time-transformed explicit symplectic integrator (TSI) or “Algorithmic regularization” (AR), which can be used for any potential that only depends on the coordinates of particles.

Here we describe the TSI (AR) method in a bit more detail. The symplectic integrator can conserve the Hamiltonian and the angular momentum of a system. Thus it is very suitable for simulating the long-term evolution of a system. However, it requires a constant integration step. In the classical symplectic method, time step, $\mathrm{d}t$ , is also the integration step, $\mathrm{d}s$ . This leads to a low efficiency in integrating an eccentric Kepler orbit. In order to be accurately enough, $\mathrm{d}t$ has to be fixed to the smallest value determined at the pericenter. The solution is to apply a time transformation, $\mathrm{d}t=g\mathrm{d}s$ , which decouples $\mathrm{d}s$ and $\mathrm{d}t$ with a function $g$ (e.g Hairer, 1997). Thus, $\mathrm{d}t$ can vary to avoid the issue of low efficiency while $\mathrm{d}s$ is fixed to keep the symplectic property. This method can be described by the extended phase space Hamiltonian, where $t$ is treated as a coordinate and need to be integrated. The disadvantage is that the time transformation usually results in an inseparable Hamiltonian and only the expensive implicit integrator can be used. Mikkola & Tanikawa (1999) and Preto & Tremaine (1999) find a solution by defining a specific type of $g$ (see Section 3) that the Hamiltonian can be written in a separable style, in order to use the explicit symplectic integrator. Especially, for an isolated binary system, if $\mathrm{d}t=\mathrm{d}s/|U|$ , where $|U|$ is the absolute value of the binary potential energy (similar like the Burdet-Heggie time transformation), the integrator behaviors dramatically well for the Kepler orbit, i.e., the numerical trajectory follows the exact one with a phase error of time.

In this work, we develop an efficient method for simulating few-body systems by combining the slow-down and the TSI (AR) schemes. In Section 2, we show the slow-down Hamiltonian and the equation of motion for several types of few-body systems, such as perturbed binaries, triples and general hierarchical systems. We demonstrate that the slow-down method can correctly reproduce the secular evolution and conserve the angular momentum. In Section 3, the TSI (AR) method is described in detail. The combination of two, the slow-down time-transformed symplectic method, is discussed in Section 4. The implementation and numerical tests are shown in Section 5 and 6. Finally, we discuss our results and draw conclusions in Section 7.

To be convenient, hereafter $\bm{r}$ , $\bm{p}$ , $\bm{v}$ and $m$ represent coordinates, conjugate momenta, velocities and masses of particles separately. We define the phase-space vector as $\bm{w}\equiv(\bm{r},\bm{p})$ . The suffix of a variable using an index (e.g. $i$ , $j$ , $1$ , $2$ …) represents a specific particle while no suffix represents all particles in a $N$ -body system. For exmample, $\bm{w}\equiv$ ( $\bm{w}_{1}$ , $\bm{w}_{2}$ ,…, $\bm{w}_{N}$ ).

2 Slow-down Hamiltonian

The slow-down method can be described by a modified Hamiltonian (Mikkola & Aarseth, 1996). For a perturbed binary, the slow-down Hamiltonian is

H_{\mathrm{sd}}=\frac{1}{\kappa}H_{\mathrm{b}}+(H-H_{\mathrm{b}}),

(1)

where $H$ is the original Hamiltonian of the system, $H_{\mathrm{b}}$ is the Hamiltonian of the binary and $\kappa$ is the slow-down factor. The equation of motion is

\displaystyle\frac{\mathrm{d}\bm{w}}{\mathrm{d}t}=\{\bm{w},H_{\mathrm{sd}}\},

(2)

where $\{\}$ is the Possion bracket. When $\kappa=1$ , $H_{\mathrm{sd}}\equiv H$ . The important point is that although the Hamiltonian is modified, $\bm{w}$ keeps the original definition. Thus, the slow-down affects the time evolution of orbital phase, but the positions and velocities of particles are not scaled.

2.1 Perturbed binary system

We first consider a simple example of one binary with an external potential $U(\bm{r})$ and a fixed $\kappa$ ,

H_{\mathrm{sd}}=\frac{1}{\kappa}\left[\frac{1}{2}m_{1}\bm{v}_{1}^{2}+\frac{1}{2}m_{2}\bm{v}_{2}^{2}-\frac{Gm_{1}m_{2}}{|\bm{r}_{1}-\bm{r}_{2}|}\right]+U(\bm{r}).

(3)

where the suffixes “1” and “2” denote the two components of the binary, and $\bm{p}_{i}$ is replaced by $m_{i}\bm{v}_{i}$ .

Applying Eq. 2, the evolution of $\bm{r}$ and $\bm{v}$ can be described by

	$\displaystyle\frac{\mathrm{d}\bm{v}_{i}}{\mathrm{d}t}=$	$\displaystyle-\frac{1}{\kappa}\frac{Gm_{3-i}(\bm{r}_{i}-\bm{r}_{3-i})}{\|\bm{r}_{i}-\bm{r}_{3-i}\|^{3}}-\frac{\partial U(\bm{r})}{\partial\bm{r}_{i}}$	$\displaystyle(i=1~{}\mathrm{or}~{}2)$		(4)
	$\displaystyle\frac{\mathrm{d}\bm{r}_{i}}{\mathrm{d}t}=$	$\displaystyle\frac{1}{\kappa}\bm{v}_{i}.$			(4)

Thus, in the reference frame of the perturber, the motion of binary is slowed down by a factor of $\kappa$ , while in the view of the binary, the perturbation is enhanced by $\kappa$ times. We can understand this better by considering the cumulative effect of perturbation to $\bm{v}_{i}$ on one slow-down orbit.

The effective period of the slow-down binary has $P_{\mathrm{sd}}=\kappa P$ , where $P$ is the original period. By integrating $\mathrm{d}\bm{v}_{i}/\mathrm{d}t$ of Eq. 4, the change of $\bm{v}_{i}$ after one $P_{\mathrm{sd}}$ is

\delta\bm{v}_{i}=\int_{0}^{P_{\mathrm{sd}}}{\frac{\mathrm{d}\bm{v}_{i}}{\mathrm{d}t}}.

(5)

In the case of no perturbation, $\bm{v}_{i}(t=P_{\mathrm{sd}})=\bm{v}_{i}(t=0)$ . Thus

\int_{0}^{P_{\mathrm{sd}}}{\frac{\mathrm{d}\bm{v}_{i}}{\mathrm{d}t}}=0.

(6)

When the effect of $U(\bm{r})$ is weak,

\delta\bm{v}_{i}\approx\int_{0}^{P_{\mathrm{sd}}}{-\frac{\partial U(\bm{r})}{\partial\bm{r}_{i}}}.

(7)

Then $\delta\bm{v}_{i}$ is only determined by the integrated effect of perturbation. Here $\bm{r}_{i}$ only passes one cycle, but the integration time is $\kappa P$ . Thus, the effect the perturbation is amplified by $\kappa$ times on one orbit of binary.

On the other hand, when $\kappa$ is an integer, we can measure $\delta\bm{v}_{i}$ in the original case (no slow-down; denoted as $\delta\bm{v^{\prime}}_{i}$ ) after the same physical time by a summation of $\kappa$ pieces of integral, where each piece represents one orbit:

	$\displaystyle\delta\bm{v^{\prime}}_{i}$	$\displaystyle=\sum_{j=1}^{\kappa}\int_{(j-1)P}^{jP}{-\frac{\partial U(\bm{r^{\prime}})}{\partial\bm{r^{\prime}}_{i}}}$		(8)
		$\displaystyle=\sum_{j=1}^{\kappa}\mathcal{V}_{j},$		(8)

If the perturbation does not change significantly,

\displaystyle\mathcal{V}_{j}\approx\mathcal{V}_{1}

\displaystyle(j=1,\kappa).

(9)

This means the effect of perturbation can be represented by the integration of one orbit ( $\mathcal{V}_{1}$ ) multiplied by $\kappa$ times. Thus, it is equivalent to the treatment of the slow-down method and $\delta\bm{v}_{i}\approx\delta\bm{v^{\prime}}_{i}$ . In other words, the slow-down method approximates the perturbation in an orbit averaged way. However, when the external potential changes significantly within the time interval of $\kappa P$ , Eq. 9 becomes invalid. Thus, there is a threshold of $\kappa$ . We discuss the criterion in Section 2.4.

2.2 Triple system

For a hierarchical triple,

	$\displaystyle H_{\mathrm{sd}}=$	$\displaystyle\frac{1}{\kappa}\left[\frac{1}{2}m_{1}(\bm{v}_{1}-\bm{v}_{\mathrm{cm}})^{2}+\frac{1}{2}m_{2}(\bm{v}_{2}-\bm{v}_{\mathrm{cm}})^{2}-\frac{Gm_{1}m_{2}}{\|\bm{r}_{1}-\bm{r}_{2}\|}\right]$		(10)
	$\displaystyle+$	$\displaystyle\frac{1}{2}(m_{1}+m_{2})\bm{v}_{\mathrm{cm}}^{2}+\frac{1}{2}m_{3}\bm{v}_{3}^{2}-\frac{Gm_{1}m_{3}}{\|\bm{r}_{1}-\bm{r}_{3}\|}-\frac{Gm_{2}m_{3}}{\|\bm{r}_{2}-\bm{r}_{3}\|},$		(10)

where the suffixes “1” and “2” denote the two components of the inner binary and “3” represents the third body, and $\bm{v}_{\mathrm{cm}}$ is the center-of-the-mass velocity of the inner binary. Because the slow-down affects only the internal motion of the binary in its rest frame, $\bm{v}_{\mathrm{cm}}$ is subtracted from the slow-down term.

Applying Eq. 2 (with a fixed $\kappa$ to Eq. 10), the evolution of $\bm{r}$ and $\bm{v}$ can be described by

•

binary components ( $i=1~{}\mathrm{or}~{}2$ ):

	$\displaystyle\frac{\mathrm{d}\bm{v}_{i}}{\mathrm{d}t}=$	$\displaystyle-\frac{1}{\kappa}\frac{Gm_{3-i}(\bm{r}_{i}-\bm{r}_{3-i})}{\|\bm{r}_{i}-\bm{r}_{3-i}\|^{3}}-\frac{Gm_{3}(\bm{r}_{i}-\bm{r}_{3})}{\|\bm{r}_{i}-\bm{r}_{3}\|^{3}}$		(11)
	$\displaystyle\frac{\mathrm{d}\bm{r}_{i}}{\mathrm{d}t}=$	$\displaystyle~{}\frac{1}{\kappa}(\bm{v}_{i}-\bm{v}_{\mathrm{cm}})+\bm{v}_{\mathrm{cm}};$		(11)

•

third body ( $i=3$ ):

	$\displaystyle\frac{\mathrm{d}\bm{v}_{i}}{\mathrm{d}t}=$	$\displaystyle-\sum_{j=1}^{2}\frac{Gm_{j}(\bm{r}_{i}-\bm{r}_{j})}{\|\bm{r}_{i}-\bm{r}_{j}\|^{3}}$		(12)
	$\displaystyle\frac{\mathrm{d}\bm{r}_{i}}{\mathrm{d}t}=$	$\displaystyle~{}\bm{v}_{i}.$		(12)

The form of the third body is identical to the original one.

Refer to caption — Figure 1: The illustration to show in a hierarchical triple how the force from the perturber to the binary component is calculated in the slow-down and the original methods. The left circle represents the orbit of one binary component and the right curve indicates the track of the perturber. The points alone the two tracks indicate the corresponding positions (same color or number) for calculating the perturbation force. The number represents the original method where binary pass two cycles while the color represents the slow-dowm method with $\kappa=2$ (one cycle). For example, the two solid lines and one dashed line link the corresponding positions (dark blue point with number of $7$ ) in the original and the slow-down cases separately.

We demonstrate how the perturbation force is calculated in the slow-down and the original integration methods in Fig. 1. Along the binary orbit (one of the component), $8$ sample points with an equal angular separation are chosen for the demonstration. When the binary finish two cycles in the original case, each sample point has twice force calculations from the perturber at the corresponding positions shown with the same number along its track. In the slow-down case with $\kappa=2$ , the binary finish one orbit when the perturber passes the same track. Thus, the corresponding position of the perturber (with same colors) for each sample point is different from that of the original case. If the sample points indicate the step sizes in an integrator, the slow-down method reduces the step sizes by $\kappa$ times.

The two types of lines that connect the dark blue point $7$ represent one example of the perturbation calculation in the two methods separately. If the solid and dashed lines have very close length and direction, the slow-down method give a good approximation of the perturbation force. This is the case when the perturber is far away ( $|\bm{r}_{3}|\gg|\bm{r}_{1}-\bm{r}_{2}|$ ), i.e., the perturbation is weak.

2.2.1 Orbit averaged Hamiltonian

For a hierarchical triple, we can also describe the slow-down method in terms of the orbit averaged method. By using the Legendre expansion for the perturbation term (e.g. Harrington, 1968; Naoz, et al., 2013; Naoz, 2016), $H_{\mathrm{sd}}$ can be written as the combination of three components:

$\displaystyle H_{\mathrm{sd}}=$	$\displaystyle\frac{1}{\kappa}H_{\mathrm{b,in}}+H_{\mathrm{b,out}}+H_{\mathrm{pert}}$	(13)
$\displaystyle H_{\mathrm{b,in}}=$	$\displaystyle-\frac{Gm_{1}m_{2}}{2a_{\mathrm{in}}}$
$\displaystyle H_{\mathrm{b,out}}=$	$\displaystyle-\frac{G(m_{1}+m_{2})m_{3}}{2a_{\mathrm{out}}}$
$\displaystyle H_{\mathrm{pert}}=$	$\displaystyle-\frac{G}{\|\bm{r}_{\mathrm{out}}\|}\sum_{n=2}^{\infty}\mathcal{M}_{n}\left(\frac{\|\bm{r}_{\mathrm{in}}\|}{\|\bm{r}_{\mathrm{out}}\|}\right)^{n}P_{n}(\cos{\Psi}),$

where $P_{n}$ is the Legendre polynomials, the suffixes “in” and “out” indicate the inner and outer binaries, $a$ represents the semi-major axis, $\bm{r}_{\mathrm{in}}=\bm{r}_{2}-\bm{r}_{1}$ , $\bm{r}_{\mathrm{out}}=\bm{r}_{3}-\bm{r}_{\mathrm{cm}}$ , $\bm{r}_{\mathrm{cm}}$ is the center-of-the-mass position of the inner binary and $\Psi$ is the angle between $\bm{r}_{\mathrm{in}}$ and $\bm{r}_{\mathrm{out}}$ . The mass coefficient sequence,

\mathcal{M}_{n}=m_{1}m_{2}m_{3}\frac{m_{1}^{n-1}-(-m_{2})^{n-1}}{(m_{1}+m_{2})^{n}}.

(14)

Using the Delaunay’s elements, the two binary terms of Hamiltonian have the form

	$\displaystyle H_{\mathrm{b,in}}=$	$\displaystyle-\frac{G^{2}m_{1}^{3}m_{2}^{3}}{2(m_{1}+m_{2})\mathcal{L}_{\mathrm{in}}^{2}}$		(15)
	$\displaystyle H_{\mathrm{b,out}}=$	$\displaystyle-\frac{G^{2}(m_{1}+m_{2})^{3}m_{3}^{3}}{2(m_{1}+m_{2}+m_{3})\mathcal{L}_{\mathrm{out}}^{2}}$		(15)

where the conjugate momenta,

	$\displaystyle\mathcal{L}_{\mathrm{in}}=$	$\displaystyle\frac{m_{1}m_{2}}{m_{1}+m_{2}}\sqrt{G(m_{1}+m_{2})a_{\mathrm{in}}}$		(16)
	$\displaystyle\mathcal{L}_{\mathrm{out}}=$	$\displaystyle\frac{(m_{1}+m_{2})m_{3}}{m_{1}+m_{2}+m_{3}}\sqrt{G(m_{1}+m_{2}+m_{3})a_{\mathrm{out}}}.$		(16)

Applying the equation of motion, the time derivations of their corresponding coordinates, $\ell_{\mathrm{in}}$ and $\ell_{\mathrm{out}}$ , are

	$\displaystyle\frac{\mathrm{d}\ell_{\mathrm{in}}}{\mathrm{d}t}=$	$\displaystyle\frac{\partial H_{\mathrm{sd}}}{\partial\mathcal{L}_{\mathrm{in}}}=\frac{1}{\kappa}\frac{G^{2}m_{1}^{3}m_{2}^{3}}{(m_{1}+m_{2})\mathcal{L}_{\mathrm{in}}^{3}}+\frac{\partial H_{\mathrm{pert}}}{\partial\mathcal{L}_{\mathrm{in}}}$		(17)
	$\displaystyle\frac{\mathrm{d}\ell_{\mathrm{out}}}{\mathrm{d}t}=$	$\displaystyle\frac{\partial H_{\mathrm{sd}}}{\partial\mathcal{L}_{\mathrm{out}}}=\frac{G^{2}(m_{1}+m_{2})^{3}m_{3}^{3}}{(m_{1}+m_{2}+m_{3})\mathcal{L}_{\mathrm{out}}^{3}}+\frac{\partial H_{\mathrm{pert}}}{\partial\mathcal{L}_{\mathrm{out}}}.$		(17)

The effect of slow-down is reflected on the first term of $\mathrm{d}\ell_{\mathrm{in}}/\mathrm{d}t$ , which represents the mean motion of the inner binary. Without the perturbation term, Eq. 17 indicates that the mean motion of the inner binary is $\kappa$ times slower of that in the original case, as discussed in Section 2.1. The orbit averaged method uses a canonical transformation (Von Zeipel transformation) to remove the short-period terms ( $\ell_{\mathrm{in}}$ and $\ell_{\mathrm{out}}$ ) in $H_{\mathrm{pert}}$ (e.g. Naoz, et al., 2013):

\overline{H_{\mathrm{pert}}}=\int_{0}^{2\pi}{\int_{0}^{2\pi}{H_{\mathrm{pert}}\mathrm{d}\ell_{\mathrm{in}}\mathrm{d}\ell_{\mathrm{out}}}}.

(18)

Since $\ell_{\mathrm{in}}$ and $\ell_{\mathrm{out}}$ do not appear in $H_{\mathrm{b,in}}$ and $H_{\mathrm{b,out}}$ and $H_{\mathrm{pert}}$ is independent of (the fixed) $\kappa$ , the transformation results in the same form of $\overline{H_{\mathrm{pert}}}$ in the slow-down and original cases. Thus, the slow-down method provides the same secular evolution as the orbit averaged way.

2.3 $N$ -body system with multiple slow-down binaries

For a general few-body system containing $N_{\mathrm{b}}$ binaries, each binary can have an individual slow-down factor, $\kappa_{i}$ . The slow-down Hamiltonian has the general form:

$\displaystyle H_{\mathrm{sd}}=$	$\displaystyle\sum_{i}^{N_{\mathrm{b}}}\frac{1}{\kappa_{i}}H_{\mathrm{b},i}+\left(H-\sum_{i}H_{\mathrm{b},i}\right)$	(19)
$\displaystyle H_{\mathrm{b},i}=$	$\displaystyle\frac{1}{2}m_{1,i}(\bm{v}_{1,i}-\bm{v}_{\mathrm{cm},i})^{2}+\frac{1}{2}m_{2,i}(\bm{v}_{2,i}-\bm{v}_{\mathrm{cm},i})^{2}$
	$\displaystyle-\frac{Gm_{1,i}m_{2,i}}{\|\bm{r}_{1,i}-\bm{r}_{2,i}\|}$
$\displaystyle\bm{v}_{\mathrm{cm},i}=$	$\displaystyle\frac{m_{1,i}\bm{v}_{1,i}+m_{2,i}\bm{v}_{2,i}}{m_{1,i}+m_{2,i}},$

where the suffixes “ $1,i$ ” and “ $2,i$ ” represent the two components in the $i^{\mathrm{th}}$ binary. To be convenient, we use the first $2N_{\mathrm{b}}$ indices to indicate binary components and others to indicate singles, i.e., “1,i” and “2,i” are equivalent to $2i-1$ and $2i$ separately.

Applying Eq. 2 with all $\kappa_{i}$ being fixed, the equation of motion is

•

binary components ( $i=1,2,...,N_{\mathrm{b}};k=1~{}\mathrm{or}~{}2$ ):

$\displaystyle\frac{\mathrm{d}\bm{v}_{k,i}}{\mathrm{d}t}=$	$\displaystyle-\frac{1}{\kappa_{i}}\frac{Gm_{3-k,i}(\bm{r}_{k,i}-\bm{r}_{3-k,i})}{\|\bm{r}_{k,i}-\bm{r}_{3-k,i}\|^{3}}$	(20)
	$\displaystyle-\sum_{j=1;j\neq 2i+1-k}^{N}\frac{Gm_{j}(\bm{r}_{k,i}-\bm{r}_{j})}{\|\bm{r}_{k,i}-\bm{r}_{j}\|^{3}}$
$\displaystyle\frac{\mathrm{d}\bm{r}_{k,i}}{\mathrm{d}t}=$	$\displaystyle~{}\frac{1}{\kappa_{i}}(\bm{v}_{k,i}-\bm{v}_{\mathrm{cm},i})+\bm{v}_{\mathrm{cm},i};$

•

single body ( $i=2N_{\mathrm{b}}+1,2N_{\mathrm{b}}+2,...,N$ ):

	$\displaystyle\frac{\mathrm{d}\bm{v}_{i}}{\mathrm{d}t}=$	$\displaystyle-\sum_{j=1;j\neq i}^{N}\frac{Gm_{j}(\bm{r}_{i}-\bm{r}_{j})}{\|\bm{r}_{i}-\bm{r}_{j}\|^{3}}$		(21)
	$\displaystyle\frac{\mathrm{d}\bm{r}_{i}}{\mathrm{d}t}=$	$\displaystyle~{}\bm{v}_{i}.$		(21)

This general form suggests that $\kappa_{i}$ only explicitly affects the internal motion of the binary $i$ . The force calculations of singles and other binaries do not explicitly depend on $\kappa_{i}$ .

2.4 Varying $\kappa$ in the integration

In previous sections, we assume that $\kappa$ is constant. In many few-body systems, the perturbation to a binary can vary significantly, especially when a nearby perturber has a highly eccentric or hyperbolic orbit. Thus, varying $\kappa$ during the integration is necessary to balance the performance and accuracy. Mikkola & Aarseth (1996) suggests to calculate $\kappa$ based on the ratio between the tidal force from perturbers and the internal force of the binary assuming the separation is $2a_{\mathrm{in}}$ . For a system with one binary and $N_{\mathrm{p}}$ perturbers, we define a modified version by including the eccentricity ( $e$ ),

\kappa(\bm{w})=k_{\mathrm{ref}}\frac{m_{1}m_{2}}{(m_{1}+m_{2})\left[a_{\mathrm{in}}(1+e_{\mathrm{in}})\right]^{3}}\sum_{i}^{N_{\mathrm{p}}}\frac{|\bm{r}_{i}-\bm{r}_{\mathrm{cm}}|^{3}}{m_{i}},

(22)

where $k_{\mathrm{ref}}$ is a constant coefficient.

Since this $\kappa$ depends on $\bm{w}$ , one might think that the additional terms of ( $H_{\mathrm{b,in}}\partial\kappa(\bm{w})/\partial\bm{r}$ , $H_{\mathrm{b,in}}\partial\kappa(\bm{w})/\partial\bm{v}$ ) should appear in the equation of motion. However, these terms actually should not be included to ensure a correct secular evolution. In the case of a triple system, by applying Eq. 22 in the Hamiltonian form of Eq. 13,

\frac{1}{\kappa}H_{\mathrm{b,in}}=-\frac{k_{\mathrm{ref}}G(m_{1}+m_{2})m_{3}}{2|\bm{r}_{\mathrm{out}}|}\left[\frac{a_{\mathrm{in}}(1+e_{\mathrm{in}})}{|\bm{r}_{\mathrm{out}}|}\right]^{2}.

(23)

$H_{\mathrm{b,in}}$ becomes to depend on $\ell_{\mathrm{out}}$ (from $a_{\mathrm{in}}/|\bm{r}_{\mathrm{out}}|$ ) in the same order as that of the perturbation term in $H_{\mathrm{pert}}$ with $n=2$ . The additional terms including $H_{\mathrm{b,in}}$ then affects the orbital average of $\ell_{\mathrm{out}}$ (Eq. 18), which may not result in the same secular evolution as the orbit averaged method.

On the other hand, when Eq. 9 is satisfied, $\kappa(\bm{w})$ is expected to evolve slowly in one step. As the first order approximation, $\kappa$ can be treated as a temporary constant and is only updated at the end of one integration step. Essentially, $\kappa$ is only a step function of time. Then the unwarranted effect can be avoided. We can describe this Jumping- $\kappa$ method for a system with $N_{\mathrm{b}}$ binaries as:

1.

advance $\mathrm{d}t$ with fixed $\kappa_{i}$ measured at $t$ .
2.

update $\kappa_{i}$ at the new time $t^{\prime}=t+\mathrm{d}t$ ;

apply an instant correction of $H_{\mathrm{sd}}$ by

\delta H_{\mathrm{sd}}(t^{\prime})=\sum_{i}^{N_{\mathrm{b}}}{\left(\frac{1}{\kappa_{i}(t^{\prime})}-\frac{1}{\kappa_{i}(t)}\right)H_{\mathrm{b},i}(t^{\prime})}.

(24)

With the step (iii), even $H_{\mathrm{sd}}$ is not conserved, the cumulative numerical error of $H_{\mathrm{sd}}$ can still be properly tracked as

\epsilon(H_{\mathrm{sd}},k\mathrm{d}t)=H_{\mathrm{sd}}(k\mathrm{d}t)-\sum_{i=1}^{k}\delta H_{\mathrm{sd}}(i\mathrm{d}t)-H_{\mathrm{sd}}(0),

(25)

where $k$ is the total step count.

In this method, although the unwarranted effect on the secular evolution is not included in one integration step, the instant change of $H_{\mathrm{sd}}$ may still introduce an additional effect. Thus, it is necessary to limit the change rate of $\kappa$ for safety. We estimate how $\kappa(\bm{w})$ changes depending on the orbital parameters. With Eq. 22, the differential of $1/\kappa$ ,

	$\displaystyle\mathrm{d}\left(\frac{1}{\kappa}\right)=$	$\displaystyle k_{\mathrm{ref}}\sum_{i}^{N_{\mathrm{p}}}\frac{3(m_{1}+m_{2})m_{i}}{m_{1}m_{2}}\frac{a^{3}}{\|\bm{r}_{i}-\bm{r}_{\mathrm{cm}}\|^{3}}$		(26)
		$\displaystyle\left(\frac{\mathrm{d}a}{a}-\frac{\mathrm{d}\|\bm{r}_{i}-\bm{r}_{\mathrm{cm}}\|}{\|\bm{r}_{i}-\bm{r}_{\mathrm{cm}}\|}\right).$		(26)

Typically, in the weak perturbation condition, $a$ does not change significantly but $|\bm{r}_{i}-\bm{r}_{\mathrm{cm}}|$ may vary a lot if the perturber has an eccentric or hyperbolic orbit. In the case of one perturber, Eq. 26 can be simplified as

\mathrm{d}\left(\frac{1}{\kappa}\right)\approx-3k_{\mathrm{ref}}\left(\frac{1}{\kappa}\right)\frac{\mathrm{d}|\bm{r}_{3}-\bm{r}_{\mathrm{cm}}|}{|\bm{r}_{3}-\bm{r}_{\mathrm{cm}}|}.

(27)

Thus, the time derivation of $1/\kappa$ depends on $|\bm{v}_{3}-\bm{v}_{\mathrm{cm}}|/|\bm{r}_{3}-\bm{r}_{\mathrm{cm}}|$ . The evolution timescale is the order of the crossing time of the perturber.

The slow change of $\kappa$ requires that the strength of perturbation does not vary significantly within at least one $P_{\mathrm{sd}}$ . Based on Eq. 27, we can introduce a timescale criterion to limit the maximum value of $\kappa$ for safety:

\kappa_{\mathrm{max}}=\frac{c|\langle\bm{r}_{i}-\bm{r}_{3}\rangle|}{P_{\mathrm{sd}}|\langle\bm{v}_{i}-\bm{v}_{\mathrm{cm}}\rangle|},

(28)

where $c$ is a constant coefficient; $\langle\bm{r}_{i}-\bm{r}_{3}\rangle$ and $\langle\bm{v}_{i}-\bm{v}_{\mathrm{cm}}\rangle$ represent the mass weighted average of velocities and positions of all perturbers referring to the center-of-the-mass of the binary separately. If a binary is the perturber, the mass weights help to remove the velocity fluctuations caused by the internal motion of the two components of the perturber. For the Jumping- $\kappa$ method, Eq. 22 and 28 together provide the criterion to estimate $\kappa$ . In Section 6, we show that the numerical tests by using the Jumping- $\kappa$ method indeed provide a correct secular evolution.

2.5 Conserved quantities

2.5.1 Energy

It is important to know, with the irregular form of $H_{\mathrm{sd}}$ , what are conserved quantities during the evolution that can be used to check the quality of the integration. When $\kappa$ is a constant, $H_{\mathrm{sd}}$ does not explicitly depend on time. Thus the “slow-down energy”, $H_{\mathrm{sd}}$ , is a conserved quantity, which means that the physical energy, $H$ , changes during the evolution. The variation of $H$ can be calculated by

\delta H=H-H_{\mathrm{sd}}=\sum_{i}^{N_{\mathrm{b}}}{\left(1-\frac{1}{\kappa_{i}}\right)H_{\mathrm{b},i}}.

(29)

Thus, $H$ oscillates in a timescale of the secular evolution of the binaries.

When $\kappa$ is not fixed, whether $H_{\mathrm{sd}}$ is conserved depends on how $\kappa$ is evaluated. If $\kappa$ follows Eq. 22 exactly and the equation of motion contains the additional terms of ( $H_{\mathrm{b,in}}\partial\kappa(\bm{w})/\partial\bm{r}$ , $H_{\mathrm{b,in}}\partial\kappa(\bm{w})/\partial\bm{v}$ ), $\kappa$ would not explicitly depend on time, thus $H_{\mathrm{sd}}$ is conserved. However, this also indicates that a significant variation of $\kappa$ results in a large change of binary energy. Such behaviour is introduced by the unwarranted effect discussed in Section 2.4. Instead, the Jumping- $\kappa$ method that exclude the additional term can avoid such problem, although $H_{\mathrm{sd}}$ explicitly depends on time and is not conserved. Since the physical energy is anyway not conserved, it is more important to ensure the correct behaviour of secular evolution rather than to conserve $H_{\mathrm{sd}}$ . On the other hand, we can still follow the integration error by using Eq. 25 in the absence of the energy conservation.

2.5.2 Angular momentum

The angular momentum defined by

\bm{L}=\sum_{i}m_{i}\bm{r}_{i}\times\bm{v}_{i}

(30)

is a conserved quantity in the original case since

\{\bm{L},H\}=\bm{0},

(31)

where $\bm{0}$ is a zero vector. Interestingly, it can be proved that $\bm{L}$ is also conserved with $H_{\mathrm{sd}}$ (Eq. 19; see Appendix), i.e.

\{\bm{L},H_{\mathrm{sd}}\}=\bm{0}.

(32)

This conservation is not only valid for an invariant $\kappa$ . With the help of a symbolic calculator ¹¹1We use python.sympy to prove Eq. 32., it can be shown that for $\kappa$ defined in Eq. 22, Eq. 32 is also satisfied. On the other hand, if $\kappa$ explicitly depends on time, $\bm{L}$ is still conserved because the differential in Eq. 32 is independent of time. Thus, $\bm{L}$ is generally conserved in the slow-down method and can be used to validate the quality of the integration.

3 Time transformed symplectic (AR) method

The TSI (AR) method is accurate for solving the long-term evolution of few-body systems. It can be described by the extended phase-space Hamiltonian (e.g. Hairer, 1997; Preto & Tremaine, 1999):

\Gamma(\bm{W})=g(\bm{W})[H(\bm{w},t)-H(\bm{w}(0),0)],

(33)

where $H(\bm{w},t)$ is the standard Hamiltonian and $g(\bm{W})$ is time-transformation function. The extended phase-space vector, $\bm{W}=(\bm{R},\bm{P})$ , is defined by including an additional pair of the coordinate, $t$ , and the conjugate momentum, $p_{\mathrm{t}}=-H(\bm{w}(0),0)$ (negative initial energy). $\bm{w}(0)$ is the initial value of $\bm{w}$ .

With the new differential variable $s$ ,

g(\bm{W})=\frac{\mathrm{d}t}{\mathrm{d}s}.

(34)

A special type of $g(\bm{W})$ introduced by Preto & Tremaine (1999),

g(\bm{W})=\frac{f(T(\bm{P}))-f(-U(\bm{R}))}{T(\bm{P})+U(\bm{R})},

(35)

leads to a separable $\Gamma$ :

\Gamma(\bm{W})=f(T(\bm{P}))-f(-U(\bm{R})),

(36)

where the kinetic energy, $T(\bm{P})\equiv T(\bm{p})+p_{\mathrm{t}}$ , and the potential energy, $U(\bm{R})\equiv U(\bm{r},t)$ .

Putting $\Gamma(\bm{W})$ in the equation of motion,

\frac{\mathrm{d}\bm{W}}{\mathrm{d}s}=\{\bm{W},\Gamma(\bm{W})\},

(37)

the derivatives of $\bm{W}$ with respect to $s$ are

\begin{split}\frac{\mathrm{d}\bm{r}_{i}}{\mathrm{d}s}&=f^{\prime}(T(\bm{p})+p_{\mathrm{t}})\frac{\partial T(\bm{p})}{\partial\bm{p}_{i}}\\ \frac{\mathrm{d}t}{\mathrm{d}s}&=f^{\prime}(T(\bm{p})+p_{\mathrm{t}})\\ \frac{\mathrm{d}\bm{p}_{i}}{\mathrm{d}s}&=f^{\prime}(-U(\bm{r},t))\frac{\partial U(\bm{r},t)}{\partial\bm{r}_{i}}\\ \frac{\mathrm{d}p_{\mathrm{t}}}{\mathrm{d}s}&=f^{\prime}(-U(\bm{r},t))\frac{\partial U(\bm{r},t)}{\partial t}.\\ \end{split}

(38)

Preto & Tremaine (1999) and Mikkola & Tanikawa (1999) showed that by choosing

f(x)=\log{x},

(39)

the Leap-Frog method with drift-kick-drift (DKD) mode can ensure that the numerical trajectory of a Kepler orbit follows the exact one, $\bm{w}(t)$ , with only a phase error of time. Hereafter we use “LogH” to represent this algorithm. If the time error is not important, the LogH method is very efficient to integrate a weakly perturbed Kepler orbit. Fig. 2 shows the schedule of one DKD loop.

When $H$ does not explicitly depends on $t$ , $H(\bm{w},t)=H(\bm{w}(0),0)$ , and $T(\bm{P})+U(\bm{R})=0$ . Thus $f^{\prime}(T(\bm{P}))=f^{\prime}(-U(\bm{R}))$ . Mikkola & Aarseth (2002) found that instead of calculating $T(\bm{P})$ , one can also define a variable,

u=\int\frac{\partial U(\bm{R})}{\partial\bm{R}}\cdot\frac{\mathrm{d}\bm{R}}{\mathrm{d}t}

(40)

Then $f^{\prime}(u)=f^{\prime}(T(\bm{P}))$ . The integration schedule of this method for one DKD step is shown in Fig. 3.

3.1 Geometric feature

Via the Taylor expansion, Eq. 35 can be approximated as (Preto & Tremaine, 1999)

g(\bm{W})\approx f^{\prime}(-U(\bm{R}))

(41)

By applying Eq. 39 for a Kepler orbit system, the relation between $t$ and $s$ is

\mathrm{d}t\approx\frac{\mathrm{d}s}{-U(\bm{R})}=\frac{|\bm{r}_{1}-\bm{r}_{2}|\mathrm{d}s}{Gm_{1}m_{2}}.

(42)

This is equivalent to the case in Burdet-Heggie regularization method (Burdet, 1967, 1968; Heggie, 1973). It can be shown that this time transformation can remove the singularity by including the energy and the eccentric vector in the equation of motion. This suggests that the LogH method can well handle the high-eccentric orbit.

On the other hand, we can show that $s$ has a geometric meaning. For a Kepler orbit.

|\bm{r}_{1}-\bm{r}_{2}|=a(1-e\cos E)

(43)

where $E$ is eccentric anomaly. If $E=0$ indicates time zero, the relation between $E$ and $t$ is

t=\frac{P}{2\pi}(E-e\sin E).

(44)

The differential of this equation results in

\mathrm{d}t=\frac{P}{2\pi}(1-e\cos E)\mathrm{d}E.

(45)

By using Eq. 43, the time derivative of $E$ is

\frac{\mathrm{d}E}{\mathrm{d}t}=\frac{a}{|\bm{r}_{1}-\bm{r}_{2}|}\frac{2\pi}{P}.

(46)

This suggests that $s$ represents a scaled $E$ :

\mathrm{d}s=m_{1}m_{2}\sqrt{\frac{Ga}{m_{1}+m_{2}}}\mathrm{d}E=\mathcal{L}\mathrm{d}E,

(47)

Thus, the fixed $\mathrm{d}s$ in the symplectic integrator indicates a constant step of $\mathrm{d}E$ . Since $t$ is unknown before the integration finish, it is difficult to determine $\mathrm{d}s$ to stop the integration at a given $t$ . However, the geometric meaning of $s$ provides a way to stop at a certain $E$ . This is very useful for determining the number of steps per orbit in order to control the numerical accuracy of the integration.

3.2 Conservation of energy

Although the LogH method ensures that the numerical trajectory follows the exact Kepler orbit, the numerical $H(\bm{w})$ is not perfectly conserved. Especially, when the orbit reaches the pericenter, the error of $H(\bm{w})$ is large due to the truncation of floating points. Actually, based on the definition, the exact conserved “energy” is $\Gamma(\bm{W})$ . With Eq. 33, 39 and 41,

\Gamma(\bm{W})\approx\frac{|\bm{r}_{1}-\bm{r}_{2}|}{Gm_{1}m_{2}}\left[H(\bm{w},t)-H(\bm{w}(0),0)\right].

(48)

Since $|\bm{r}_{1}-\bm{r}_{2}|$ reaches the minimum at the pericenter, the large error of $H(\bm{w})$ is cancelled out in $\Gamma(\bm{W})$ . In Section 6, we show the numerical test of such behaviour. This feature indicates that we should avoid the calculation the orbital parameters or stop the integration at the pericenter in order to avoid a large numerical error of the orbital elements.

4 Slow-Down time-transformed symplectic (SDAR) method

The slow-down method helps to solve the large timescale issue of a hierarchical system, while the TSI (AR) method can conserve $\Gamma$ and $\bm{L}$ . The combination of two (SDAR method) is expected to be an efficient algorithm to deal with the long-term evolution of few-body systems. In the TSI method, $\Gamma(\bm{W})$ can take any separable Hamiltonian, thus it is straight forward to implement $H_{\mathrm{sd}}(\bm{w},t)$ into $\Gamma(\bm{W})$ . With Eq. 19 and 36, we can obtain the time transformed slow-down Hamiltonian for a $N$ -body system with $N_{\mathrm{b}}$ binaries:

$\displaystyle\Gamma_{\mathrm{sd}}(\bm{W})=$	$\displaystyle f(T_{\mathrm{sd}}(\bm{P}))-f(-U_{\mathrm{sd}}(\bm{R}))$	(49)
$\displaystyle T_{\mathrm{sd}}(\bm{P})=$	$\displaystyle\sum_{i=1}^{N_{\mathrm{b}}}{\frac{1}{\kappa_{i}}T_{\mathrm{b},i}(\bm{v}_{1,i},\bm{v}_{2,i})}+\sum_{i=1}^{N_{\mathrm{b}}}T_{\mathrm{b,cm},i}(\bm{v}_{\mathrm{cm},i})$
	$\displaystyle+\sum_{i=2N_{\mathrm{b}}+1}^{N}T_{i}(\bm{v}_{i})+p_{\mathrm{t}}$
$\displaystyle U_{\mathrm{sd}}(\bm{P})=$	$\displaystyle\sum_{i=1}^{N_{\mathrm{b}}}{\frac{1}{\kappa_{i}}U_{\mathrm{b},i}(\bm{r}_{1,i},\bm{r}_{2,i})}+\left[U(\bm{r})-\sum_{i=i}^{N_{\mathrm{b}}}U_{\mathrm{b},i}(\bm{r}_{1,i},\bm{r}_{2,i})\right],$

where

$\displaystyle T_{\mathrm{b},i}(\bm{v}_{1,i},\bm{v}_{2,i})=$	$\displaystyle\frac{1}{2}m_{1,i}\left(\bm{v}_{1,i}-\bm{v}_{\mathrm{cm},i}\right)^{2}+\frac{1}{2}m_{2,i}\left(\bm{v}_{2,i}-\bm{v}_{\mathrm{cm},i}\right)^{2}$	(50)
$\displaystyle T_{\mathrm{b,cm},i}(\bm{p}_{\mathrm{cm},i})=$	$\displaystyle\frac{1}{2}(m_{1,i}+m_{2,i})\bm{v}_{\mathrm{cm},i}^{2}$
$\displaystyle U_{\mathrm{b},i}(\bm{r}_{1,i},\bm{r}_{2,i})=$	$\displaystyle-\frac{Gm_{1,i}m_{2,i}}{\|\bm{r}_{1,i}-\bm{r}_{2,i}\|}$
$\displaystyle U(\bm{r})=$	$\displaystyle-\sum_{i=1}^{N}\sum_{j=i+1}^{N}\frac{Gm_{i}m_{j}}{\|\bm{r}_{i}-\bm{r}_{j}\|}.$

The equation of motion of $\Gamma_{\mathrm{sd}}(\bm{W})$ is the same as Eq. 38 except that $T(\bm{P})$ and $U(\bm{P})$ are replaced by $T_{\mathrm{sd}}(\bm{P})$ and $U_{\mathrm{sd}}(\bm{P})$ separately.

Notice that here the conserved “energy” is $\Gamma_{\mathrm{sd}}$ (with fixed $\kappa$ ). Thus when the Jumping- $\kappa$ method (see Section 2.4) is applied, the correction term in the third step should be modified as

	$\displaystyle\delta\Gamma_{\mathrm{sd}}(t^{\prime})=$	$\displaystyle f\left[T_{\mathrm{sd}}[\bm{P}(t^{\prime}),\kappa(t^{\prime})]\right]-f\left[U_{\mathrm{sd}}[\bm{R}(t^{\prime}),\kappa(t^{\prime})]\right]-$		(51)
		$\displaystyle\left[f\left[T_{\mathrm{sd}}[\bm{P}(t^{\prime}),\kappa(t)]\right]-f\left[U_{\mathrm{sd}}[\bm{R}(t^{\prime}),\kappa(t)]\right]\right]$		(51)

The corresponding cumulative numerical error, $\varepsilon(\Gamma_{\mathrm{sd}})$ , can be obtained by replacing $H_{\mathrm{sd}}$ to $\Gamma_{\mathrm{sd}}$ in Eq. 25.

5 Implementation

Based on the SDAR method, we develop a software library, sdar, written in the C++ language with the object oriented programming ²²2URL: https://github.com/lwang-astro/SDAR. This library contains three modules: ar, hermite and kepler.

ar is the implementation of the SDAR method shown in Eq. 50 with the Jumping- $\kappa$ algorithm. The high-order explicit symplectic integrator introduced by Yoshida (1990) is used. Both the LogH and the Mikkola & Aarseth (2002) methods are implemented.

hermite is a hybrid integrator which combines a $4^{th}$ order block-time-step Hermite method and the ar module. The Hermite integrator is for the global $N$ -body system and ar deals with the compact subgroups (few-body systems). The subgroup can contain inner binaries with individual $\kappa_{i}$ , while the system as a whole can also apply the slow-down method if it has a periodical evolution. Thus, hermite allows a nested two-level slow-down treatment. This can accelerate the integration of a weakly perturbed stable hierarchical system.

kepler is a tool to construct a Hierarchical (Kepler) binary tree for a group of particles. It is used to identify the inner binaries and calculate $\kappa$ .

The quadruple-double precision library, qd (Hida, Li & Bailey, 2001), is included in the code. Thus the user can switch on the high-precision support (up to $62$ -digit precision). Notice that when qd is used, the performance is reduced and also it is not thread-safe.

6 Numerical test

In this section, we check whether the slow-down method can correctly reproduce the secular evolution of few-body systems. The comparison between the slow-down and original algorithms is done by using the ar module with the $6^{th}$ -order symplectic method (Solution A in Table 1 from Yoshida, 1990). We use the Delaunay’s elements to describe the geometric information of the binary orbits. By selecting the Cartesian coordinate system ( $x$ - $y$ - $z$ ), the three angles are (see Fig. 4):

•

$\theta$ : inclination; the angle from $z$ -axis to $\bm{L}$ .
•

$\phi$ : longitude of the ascending node (often denoted as $\Omega$ ); the angle from $x$ -axis to the intersection of the orbital plane and the $x$ - $y$ plane ( $x^{\prime}$ -axis).
•

$\psi$ : argument of periapsis (often denoted as $\omega$ ); the angle from eccentric vector to the $x^{\prime}$ -axis.

6.1 Hierarchical triple (B-S)

The Kozai-Lidov effect is one of the most important feature appearing in a family of hierarchical triple systems (Kozai, 1962; Lidov, 1962). For example, when the outer binary has a circular orbit and the secondary of the inner binary has a negligible mass compared to the other two bodies, the $z^{\prime}$ component of the orbit averaged angular momentum of the inner binary is a conserved quantity:

L_{\mathrm{z^{\prime}},2}\propto\sqrt{(1-e_{\mathrm{in}}^{2})}\cos(\theta_{\mathrm{in}})=\mathrm{const},

(52)

Where the $z$ axis is defined in the invariable plane. When $\theta_{\mathrm{in}}$ is in the range of $\cos^{-1}(\pm\sqrt{3/5})$ , significant oscillations of $e_{\mathrm{in}}$ and $\theta_{\mathrm{in}}$ can happen. In the astrophysical environment, such oscillation can result in a high eccentricity that may trigger the merger of the inner binary. Thus, it is important to validate that the slow-down method can correctly reproduce the Kozai-Lidov effect.

We perform a simulation of a hierarchical triple (B-S). The formula from Antognini (2015) is used to estimate the Kozai-Lidov oscillation timescale:

t_{\mathrm{KL}}\approx\frac{8}{15\pi}\left(1+\frac{m_{1}}{m_{3}}\right)\frac{P_{\mathrm{o}ut}^{2}}{P_{\mathrm{i}n}}(1-e_{\mathrm{out}}^{2})^{3/2}.

(53)

For the purpose of test, parameters which allow a large $\kappa$ and a short $t_{\mathrm{KL}}$ is preferred. Thus, smaller $m_{1}/m_{3}$ and $P_{\mathrm{o}ut}/P_{\mathrm{i}n}$ or higher eccentricity of outer binary are better. However, the former also reduce the maximum $\kappa$ (Eq. 22), thus the latter is good for both requirements.

One suitable initial condition is shown in Table 1, where $t_{\mathrm{KL}}\approx 45$ . For convenience, we use the scale-free unit (gravitational constant is one). This is also applied for all numerical tests discussed below.

Table 1: The initial condition of the hierarchical triple (B-S) system for test.

\mp

and

m_{\mathrm{s}}

are the masses of the primary and the secondary of inner and outer binaries. The values are shown in the scale-free unit with the gravitational constant,

G=1

	$m_{\mathrm{p}}$	$m_{\mathrm{s}}$	$a$	$e$	$\theta$	$\phi$	$\psi$	$E$	$P$
in	0.900	0.100	0.00100	0.900	1.50	0.00	0.00	3.14	$1.97\times 10^{-4}$
out	2.00	1.00	1.00	0.990	0.10	0.00	0.00	3.14	3.63

Both the inner and outer binaries are initially near the apo-centre positions ( $E\approx 3.14$ ). By choosing $k_{\mathrm{ref}}=1.0\times 10^{-6}$ as suggested by Mikkola & Aarseth (1996) in Eq. 22, $\kappa_{\mathrm{max}}\approx 52$ when the outer body is at the apo-centers. As the outer eccentricity, $e_{\mathrm{out}}$ , is large, $\kappa$ decreases to below one based on Eq. 22 at the pericenter position. This should be avoided, thus we let $\kappa\geq 1.0$ . To obtain a high numerical accuracy, we use $256$ steps per orbit of the inner binary. This results in $\Delta s\approx 6.98\times 10^{-5}$ . To reach about $4t_{\mathrm{KL}}$ , the total number of integration step without slow-down is about $2.36\times 10^{8}$ . In the slow-down case, it is about $3.8\times 10^{7}$ steps which is $6$ times smaller.

Fig. 5 shows the simulation result of the B-S system by using the LogH method with and without slow-down (hereafter named as “SD” and “ORG” methods) ³³3The scientific notation in the plotting style of $y$ -axis (described in Fig. 5) is defined by the python module matplotlib.ticker (see https://matplotlib.org/ for reference). In all figures, we follow this style in order to save space.. The Kozai-Lidov oscillation appears in the evolution of the orbital parameters. The comparison between red (SD) and black (ORG) curves clearly indicates that the SD method can well reproduce the Kozai-Lidov effect with the correct timescale. All orbital parameters, including $a$ , $e$ and three Euler angles of both inner and outer binaries are consistent with each others. The oscillation timescale is not exacted as predicted by Eq. 53. Since the secondary of the inner binary has a small but none zero mass, the B-S system does not exactly satisfy the test-particle condition that is assumed in Eq. 53.

There are a few sharp peaks in the evolution of $a$ and $e$ . This is due to the low time resolution of the plotting (only few thousands data points along the $x$ -axis). Fig. 6 shows the high-time-resolution evolution of $a_{\mathrm{in}}$ and $e_{\mathrm{in}}$ with $t$ in $0\sim 8\times 10^{-3}$ by using the ORG method. When the inner binary passes the pericenter, sharp peaks appear. The amplitude of peaks depends on the orbital phase of both inner and outer binaries. Since the time interval of one peak is very short, most of the peaks are not sampled in the low-time-resolution plot and they appear by chance. These peaks also exist in other orbital parameters but cannot be seen in the plots due to a large amplitude of the secular variation. As the peaks do not affect the secular evolution, they can be safely ignored.

$\kappa$ oscillates when the third body cycles between the apo-center and the pericenter, as shown in the upper panel of Fig. 7. The evolution of $H_{\mathrm{sd}}$ follows the pattern of $\kappa$ as described by Eq. 29. Although the scatter of $H_{\mathrm{sd}}$ is large, the behaviour of the secular evolution is correct. This indicates that whether $H_{\mathrm{sd}}$ is conserved does not matter.

However, we still can trace the numerical error by using $\varepsilon(\Gamma_{\mathrm{sd}})$ . In Fig. 8 the evolution of the cumulative numerical errors of different definitions of Hamiltonian, $\varepsilon(H)$ , $\varepsilon(H_{\mathrm{sd}})$ and $\varepsilon(\Gamma_{\mathrm{sd}})$ , are compared. $\varepsilon(H)$ is the error of original Hamiltonian, defined as $H(t)-H(0)$ . After $e_{\mathrm{in}}$ pass the maximum value for the first time, large oscillations appears on both $\varepsilon(H)$ and $\varepsilon(H_{\mathrm{sd}})$ . $\varepsilon(H)$ of the SD method has a much larger scatter compared to the ORG method, while in the case of $\varepsilon(H_{\mathrm{sd}})$ , both have a similar amplitude. This is consistent with the expectation as discussed in Section 2.5.1. On the other hand, the scatter of $\varepsilon(\Gamma)$ or $\varepsilon(\Gamma_{\mathrm{sd}})$ is two order of magnitude smaller than $\varepsilon(H_{\mathrm{sd}})$ without oscillations, which indicates that the conserved quantities are $\Gamma$ (ORG) and $\Gamma_{\mathrm{sd}}$ (SD) within one integration step. Thus, $\varepsilon(\Gamma_{\mathrm{sd}})$ represents the numerical errors of the integrator.

The angular momentum conservation is also checked in Fig. 9. The normalized cumulative error of three components,

\displaystyle\widetilde{\varepsilon}(L_{i})=\frac{L_{i}(t)-L_{i}(0)}{|\bm{L}(0)|}

\displaystyle(i=x,y,z),

(54)

and the total one,

\widetilde{\varepsilon}(\bm{L})=\frac{\sqrt{\left[\bm{L}(t)-\bm{L}(0)\right]^{2}}}{|\bm{L}(0)|},

(55)

are shown. Both the SD and the ORG methods provide a similar level of relative numerical errors ( $10^{-10}$ ) in the three components of $\bm{L}$ and also in $|\bm{L}|$ . Besides, the error is independent of the variation of $\kappa$ . This is consistent as we proved in Section 2.5.2.

6.2 Hierarchical quadruple (B-B)

The mean motion resonance cannot be properly reproduced by the slow-down method. When there are multiple binaries, the orbital resonance between inner binaries can occur if their periods are in resonance ratios. We investigate this effect by simulating a quadruple system (B-B). The initial condition (Table 2) is generated by splitting the third body of the B-S system to a binary, which has the same period as another. Here after the suffixes “in1” and “in2” indicate the two inner binaries. $t_{\mathrm{KL}}$ of the two binaries are about $45$ and $88$ . We evolve the system to $t=160$ with the same $ds$ used in the integration of the B-S system. The ORG method takes about $8.7\times 10^{8}$ steps while the SD method takes about $9.9\times 10^{7}$ steps ( $9$ times faster).

Table 2: The initial condition of the hierarchical quadruple (B-B) system.

	$m_{\mathrm{p}}$	$m_{\mathrm{s}}$	$a$	$e$	$\theta$	$\phi$	$\psi$	$E$	$P$
in1	0.900	0.100	0.00100	0.900	1.50	0.00	0.00	3.14	$1.97\times 10^{-4}$
in2	1.800	0.200	0.00126	0.900	0.10	0.00	0.00	3.14	$1.97\times 10^{-4}$
out	2.00	1.00	1.00	0.990	0.100	0.00	0.00	3.14	3.63

The orbital evolution is shown in Fig. 10. $t_{\mathrm{KL,in2}}$ of the simulated result is a much larger than the prediction of Eq. 53. The behaviour of the first inner binary is similar to the case of the B-S system. Compared to the case of B-S system, the cumulative divergence of $a$ for the SD and ORG methods becomes obvious. After about two $t_{\mathrm{KL}}$ , the relative difference is the order of $10^{-3}$ for $a_{\mathrm{in1}}$ and $a_{\mathrm{in2}}$ . However, this error can be neglected since the Kozai-Lidov effect dominates the secular evolution. Thus, the SD method can still provide a reasonable result.

On the other hand, since $a$ of two inner binaries do not evolve significantly, it is easy to observe the difference caused by the orbital resonance between the SD and the ORG method. Although the SD method cannot keep the original period ratio of the two binaries, the mean motion resonance still exist because $\kappa_{i}$ of the two binaries are related by Eq. 22. The upper panel of Fig. 11 shows the variation of $\kappa_{i}$ , where $\kappa_{\mathrm{in2}}$ is twice of $\kappa_{\mathrm{in1}}$ . Therefore, the period ratio changes from $1:1$ to $1:2$ . Since the orbits are slowed down, the resonance becomes stronger, which results in the larger drift in the SD case.

$\varepsilon(\Gamma_{\mathrm{sd}})$ keeps the order of $10^{-7}$ . Both SD and ORG methods show jumps when $e_{\mathrm{in1}}$ passes the maximum value due to the large numerical error at the pericenter of the first inner binary. Such error can be reduced with a smaller $\Delta s$ (not shown here). The behaviour of $|\bm{L}|$ is independent of $e_{\mathrm{in1}}$ and its relative error has an order of $10^{-10}$ .

6.3 Hyperbolic encounter (HB-B)

Table 3: The initial condition of the system of a hyperbolic encounter between two binaries (HB-B).

P

of the outer orbit is the time to reach the pericenter. The inner binary orbits are the same as in the B-B system except in the ORG models,

E_{\mathrm{in2}}

uses three values for testing the impact of initial orbital phases. The SD model (SD-3) chooses the third value of

E_{\mathrm{in2}}

	$m_{\mathrm{p}}$	$m_{\mathrm{s}}$	$a$	$e$	$\theta$	$\phi$	$\psi$	$E$	$P$
in1	0.900	0.100	0.00100	0.900	1.50	0.00	0.00	3.14	$1.97\times 10^{-4}$
in2	0.00900	0.00100	0.00200	0.900	1.50	0.00	0.00	3.00, 3.50, 3.14	$5.62\times 10^{-3}$
out	1.00	0.0100	-1.00	1.02500	0.100	0.00	0.00	-2.00	1.71

In star clusters, hyperbolic encounters between stars and binaries are frequent. Thus we investigate whether the slow-down method can provide an acceptable result in such case. The encounter happens in a short time interval. When a perturber is far away, the binary has a weak perturbation so that $\kappa$ is large. Once the encounter starts, $\kappa$ drops fast. Thus, it is necessary to limit the change rate of $\kappa$ , for example, by using the timescale criterion (Eq. 28).

We test a hyperbolic encounter between a massive binary and a low-mass binary (HB-B). The initial condition is shown in Table 3. The mass ratio between the two inner binaries is $100$ . The hyperbolic (outer) orbit has an initial separation of $2.86$ and a pericenter distance of $0.0250$ ( $25$ times of $a_{\mathrm{in1}}$ ). Initially $E_{\mathrm{out}}=-2.00$ so that the time to reach the pericenter is about $1.71$ . Based on Eq. 22, $\kappa_{\mathrm{in1}}$ has the maximum value of about $3\times 10^{4}$ and can reaches the minimum limit of $1.0$ . Thus, $\kappa_{\mathrm{in1}}$ varies about four order of magnitude in one encounter. By using the timescale criterion with $c=0.1$ in Eq. 28, $\kappa_{\mathrm{in1}}$ is limited to about $10^{3}$ , which is one magnitude smaller.

Due to the high mass ratio, the second inner binary feel a strong perturbation from the first inner binary. The final fate of the second binary depends on its orbital phase during the encounter. In the ORG case, three values of $E_{\mathrm{in2}}$ are used to test the effect of orbital phases. The SD model use the third value of $E_{\mathrm{in2}}$ listed in Table 3.

The result of orbital revolution is shown in Fig. 12. After the encounter, a jump appears in all orbital parameters. The initial $E_{\mathrm{in2}}$ significantly influences the final orbit of the secondary inner binary (a large divergence appears). However, the SD and ORG method (SD-3 and ORG-3) can provide a similar result when the initial $E_{\mathrm{in2}}$ are same ( $3.14$ ). Because the SD method loses the real orbital phase of the first binary, the result suggests that although the first binary is much more massive, its orbital phase has a minor impact on the evolution of the second binary. Therefore, the SD method can provide an acceptable result.

The upper panel of 13 compared $\kappa_{\mathrm{in1}}$ (for the SD-3 model) calculated by the perturbation criterion (Eq. 22) and by both the perturbation and the timescale criterion. The timescale criterion ensures that $\kappa_{\mathrm{in1}}$ changes more smoothly during the encounter. $\varepsilon(H_{\mathrm{sd}})$ has a larger error of ( $10^{-5}$ ) compared to $\varepsilon(H)$ at the beginning because $\kappa_{\mathrm{in1}}$ is large. It would be worse if the timescale criterion is not applied. Although $\varepsilon(\Gamma_{\mathrm{sd}})$ is not small, $|\bm{L}|$ is still well conserved (error has an order of $10^{-10}$ ), which again suggests that the conservation of $\bm{L}$ is independent of $\kappa$ .

With $\Delta s\approx 1.58\times 10^{-6}$ , the ORG method requires about $2.3\times 10^{8}$ steps to reach $t=4$ while the SD method only needs $3.8\times 10^{6}$ steps. Thus, the SD method provides a $60$ times faster performance. Notice here we only evolve the system in a short time interval around the time of encounter. In a star cluster, most of the time the binaries are weakly perturbed, thus averagely $\kappa_{i}\gg 10^{3}$ . Therefore, the total integration steps of binaries are significantly reduced during the long-term evolution. This example suggests that in the simulation of star clusters with many binaries, the SD method can dramatically improves the performance while the statistical result of encounters can be well reproduced.

7 Discussion and conclusion

In this work, we mathematically and numerically describe the slow-down time-transformed explicit symplectic (SDAR) method. It combines the benefit of the symplectic integrator which conserves the Hamiltonian and angular momentum and the high efficiency of the slow-down method to handle the long-term evolution of hierarchical systems and close encounters. An implementation of the method written in the c++ language, sdar, is publicly available. The code modularized for easily plugging in other $N$ -body codes.

The Hamiltonian and the corresponding equation of motion are discussed in detail. Although the physical energy ( $H$ ) is not conserved, the method can provide a correct secular evolution. We mathematically prove that $\bm{L}$ is always conserved with the slow-down method. We also discussed how to measure the numerical error, $\varepsilon(\Gamma_{\mathrm{sd}})$ , in the absence of energy conservation.

We show that in the LogH method, for a Kepler orbit the integration step $\mathrm{d}s$ has the geometric meaning of eccentric anomaly(Eq. 47). Using this feature, we can determine the number of steps per orbit by setting $ds$ using Eq. 47. This is very useful for controlling the integration error.

On the other hand, when the LogH method is implemented in a hybrid $N$ -body integrator, it is necessary to ensure that the integration can stop at a given time, in order to calculate the interactions between the members of few-body systems and the global system. However, because time is an integrated value, before the integration finishes, the next time is unknown. A few iteration of integration is needed to converge the next time to a certain value with a given limit of error. Such iteration can be expensive if the time synchronization is frequent. Thus, in order to avoid too many synchronization requests, it is necessary to design a good criterion for determination of the subsystems that need to be handled by the SDAR method. Especially, the number of iteration steps should be much less than the number of integration steps.

By using the ar code, we show three numerical examples (B-S, B-B and HB-B) to indicate that the SDAR method can well reproduce the Kozai-Lidov oscillation and can give an acceptable result of a hyperbolic encounter between two binaries. When the averaged value of $\kappa$ is large, the slow-down method can provide significant performance improvement. Especially, for binaries with weak perturbations, the method can provide a few order of magnitude less integration steps. Thus, by combining this algorithm in an hybrid integrator for simulating a large $N$ -body system, it is expected that a large fraction of binaries and hierarchical systems can be efficiently handled. Such a code will be available soon in our following-up work.

Acknowledgments

L.W. thanks the financial support from JSPS International Research Fellow (School of Science, The university of Tokyo).

References

Aarseth (2003) Aarseth S. J., 2003, Gravitational N-Body Simulations, Cambridge University Press
Antognini (2015) Antognini J. M. O., 2015, MNRAS, 452, 3610
Antonini, et al. (2016) Antonini F., Chatterjee S., Rodriguez C. L., Morscher M., Pattabiraman B., Kalogera V., Rasio F. A., 2016, ApJ, 816, 65
Askar, Szkudlarek, Gondek-Rosińska, Giersz & Bulik (2017) Askar A., Szkudlarek M., Gondek-Rosińska D., Giersz M., Bulik T., 2017, MNRAS, 464, L36
Bailyn (1995) Bailyn C. D., 1995, ARA&A, 33, 133
Banerjee, Baumgardt & Kroupa (2010) Banerjee S., Baumgardt H., Kroupa P., 2010, MNRAS, 402, 371
Breen & Heggie (2013) Breen P. G., Heggie D. C., 2013, MNRAS, 432, 2779
Binney & Tremaine (1987) Binney J., Tremaine S., 1987, Princeton, NJ, Princeton University Press
Burdet (1967) Burdet C. A., 1967, ZaMP, 18, 434
Burdet (1968) Burdet C. A., 1968, ZaMP, 19, 345
Davies, Piotto & de Angeli (2004) Davies M. B., Piotto G., de Angeli F., 2004, MNRAS, 349, 129
Di Carlo, et al. (2019) Di Carlo U. N., Giacobbo N., Mapelli M., Pasquato M., Spera M., Wang L., Haardt F., 2019, MNRAS, 487, 2947
Fujii & Portegies Zwart (2011) Fujii M. S., Portegies Zwart S., 2011, Sci, 334, 1380
Gualandris, Portegies Zwart & Sipior (2005) Gualandris A., Portegies Zwart S., Sipior M. S., 2005, MNRAS, 363, 223
Gvaramadze, Gualandris & Portegies Zwart (2009) Gvaramadze V. V., Gualandris A., Portegies Zwart S., 2009, MNRAS, 396, 570
Hairer (1997) Hairer, E., 1997, Applied Numerical Mathematics, 25, 219-227
Harrington (1968) Harrington R. S., 1968, AJ, 73, 190
Heggie (1973) Heggie D. C., 1973, ASSL, 34, ASSL…39
Heggie (1975) Heggie D. C., 1975, MNRAS, 173, 729
Heggie & Giersz (2008) Heggie D. C., Giersz M., 2008, MNRAS, 389, 1858
Hénon (1961) Hénon M., 1961, AnAp, 24, 369
Hida, Li & Bailey (2001) Hida Y., Li X. S., Bailey D. H., 2001, ”Algorithms for quad-double precision floating point arithmetic,” Proceedings 15th IEEE Symposium on Computer Arithmetic. ARITH-15 2001, Vail, CO, USA, 2001, pp. 155-162.
Hills (1975) Hills J. G., 1975, AJ, 80, 809
Hurley, et al. (2005) Hurley J. R., Pols O. R., Aarseth S. J., Tout C. A., 2005, MNRAS, 363, 293
Hypki & Giersz (2013) Hypki A., Giersz M., 2013, MNRAS, 429, 1221
Kozai (1962) Kozai Y., 1962, AJ, 67, 591
Kustaanheimo & Stiefel (1965) Kustaanheimo P., Stiefel E., 1965, J. Reine Angew. Math., 218, 204
Leigh, Sills & Knigge (2011) Leigh N., Sills A., Knigge C., 2011, MNRAS, 416, 1410
Leonard & Duncan (1988) Leonard P. J. T., Duncan M. J., 1988, AJ, 96, 222
Leonard & Duncan (1990) Leonard P. J. T., Duncan M. J., 1990, AJ, 99, 608
Lidov (1962) Lidov M. L., 1962, P&SS, 9, 719
Mikkola & Aarseth (1996) Mikkola S., Aarseth S. J., 1996, CeMDA, 64, 197
Mikkola & Tanikawa (1999) Mikkola S., Tanikawa K., 1999, MNRAS, 310, 745
Mikkola & Aarseth (2002) Mikkola S., Aarseth S., 2002, CeMDA, 84, 343
Naoz, et al. (2013) Naoz S., Farr W. M., Lithwick Y., Rasio F. A., Teyssandier J., 2013, MNRAS, 431, 2155
Naoz (2016) Naoz S., 2016, ARA&A, 54, 441
Spitzer (1987) Spitzer L., 1987, Dynamical evolution of globular clusters, Princeton University Press
O’Leary et al. (2006) O’Leary R. M., Rasio F. A., Fregeau J. M., Ivanova N., O’Shaughnessy R., 2006, ApJ, 637, 937
Portegies Zwart & McMillan (2000) Portegies Zwart S. F., McMillan S. L. W., 2000, ApJL, 528, L17
Poveda, Ruiz & Allen (1967) Poveda A., Ruiz J., Allen C., 1967, BOTT, 4, 86
Preto & Tremaine (1999) Preto M., Tremaine S., 1999, AJ, 118, 2532
Wang (2020) Wang L., 2020, MNRAS, 491, 2413
Yoshida (1990) Yoshida H., 1990, PhLA, 150, 262
Yu & Tremaine (2003) Yu Q., Tremaine S., 2003, ApJ, 599, 1129

Appendix A Proof for the conservation of angular momentum

Here we provide the proof for the conservation of angular momentum in the slow-down method. First, the conservation of angular momentum for an isolated binary can be described as

		$\displaystyle\{\bm{L}_{\mathrm{b}},H_{\mathrm{b}}\}=\bm{0}$		(56)
		$\displaystyle\bm{L}_{\mathrm{b}}=m_{1}\bm{r}_{1}\times\bm{v}_{1}+m_{2}\bm{r}_{2}\times\bm{v}_{2}$
		$\displaystyle H_{\mathrm{b}}=\frac{1}{2}m_{1}\bm{v}_{1}^{2}+\frac{1}{2}m_{2}\bm{v}_{2}^{2}-\frac{Gm_{1}m_{2}}{\|\bm{r}_{1}-\bm{r}_{2}\|}$

where $\bm{L}_{\mathrm{b}}$ , $H_{\mathrm{b}}$ are the angular momentum and Hamiltonian of a binary system in the center-of-the-mass frame.

After a slow-down factor $\kappa(t)$ is included,

\{\bm{L}_{\mathrm{b}},H_{\mathrm{sd}}\}=\{\bm{L}_{\mathrm{b}},\frac{1}{\kappa}H_{\mathrm{b}}\}=\frac{1}{\kappa}\{\bm{L}_{\mathrm{b}},H_{\mathrm{b}}\}=\bm{0}.

(57)

$\bm{L}_{\mathrm{b}}$ is still conserved. When the center-of-the-mass term is included,

	$\displaystyle H_{\mathrm{sd}}=$	$\displaystyle\frac{1}{\kappa}\left[\frac{1}{2}m_{1}(\bm{v}_{1}-\bm{v}_{\mathrm{cm}})^{2}+\frac{1}{2}m_{2}(\bm{v}_{2}-\bm{v}_{\mathrm{cm}})^{2}-\frac{Gm_{1}m_{2}}{\|\bm{r}_{1}-\bm{r}_{2}\|}\right]$		(58)
		$\displaystyle+\frac{1}{2}(m_{1}+m_{2})\bm{v}_{\mathrm{cm}}^{2},$		(58)

it can be shown that the conservation (Eq. 56) still exists.

If a new body is added to form a triple, the additional term appears in the angular momentum:

\bm{L}=\bm{L}_{\mathrm{b}}+\bm{L}_{3}.

(59)

Then

\{\bm{L},H_{\mathrm{b}}\}=\{\bm{L}_{\mathrm{b}},H_{\mathrm{b}}\}+\{\bm{L}_{3},H_{\mathrm{b}}\}=\{\bm{L}_{3},H_{\mathrm{b}}\}.

(60)

Since $\bm{L}_{3}$ is independent of binary position and velocity,

\{\bm{L},H_{\mathrm{b}}\}=\{\bm{L}_{3},H_{\mathrm{b}}\}=\bm{0}.

(61)

The additional term to $H_{\mathrm{sd}}$ from the third body is

H_{3}=\frac{1}{2}m_{3}\bm{v}_{3}^{2}-\sum_{i}^{2}{\frac{Gm_{i}m_{3}}{|\bm{r}_{i}-\bm{r}_{3}|}}.

(62)

It can be proved that

\{\bm{L},H_{3}\}=\bm{0}.

(63)

Thus,

\{\bm{L},H_{\mathrm{sd}}\}=\{\bm{L},H_{\mathrm{b}}\}+\{\bm{L},H_{3}\}=\bm{0}.

(64)

The angular momentum is conserved for the triple.

For systems with one binary and many singles, $4^{th},5^{th}...$ particles can be added one by one, and for the $k^{th}$ particle, the additional Hamiltonian term is

H_{i}=\frac{1}{2}m_{i}\bm{v}_{i}^{2}-\sum_{j}^{i-1}{\frac{Gm_{i}m_{j}}{|\bm{r}_{i}-\bm{r}_{j}|}}.

(65)

Since the second term in $H_{i}$ is a linear combination of $k-1$ pair potential energy, it is not difficult to show that

\{\bm{L},H_{i}\}=\bm{0}.

(66)

Eventually,

\{\bm{L},H_{\mathrm{sd}}\}=\bm{0}

(67)

Thus, $\bm{L}$ is conserved in such system.

On the other hand, here we add many singles to a slow-down binary and indicate that $\bm{L}$ is conserved. We can consider this process in an opposite way: add one slow-down binary to a systems of many singles, the new $\bm{L}$ is still conserved. Thus, it can be proved that adding arbitrary slow-down binaries to the systems, $\bm{L}$ is always conserved.