A slow-down time-transformed symplectic integrator for solving the few-body problem
MMon. Not. R. Astron. Soc. , 1–14 (2002) Printed 20 February 2020 (MN L A TEX style file v2.2)
A slow-down time-transformed symplectic integrator forsolving the few-body problem
Long Wang , (cid:63) , 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
Accepted – . Received –; in original form –
ABSTRACT
An accurate and efficient method dealing with the few-body dynamics is importantfor simulating collisional N -body systems like star clusters and to follow the forma-tion and evolution of compact binaries. We describe such a method which combinesthe time-transformed explicit symplectic integrator (Preto & Tremaine 1999; Mikkola& Tanikawa 1999) and the slow-down method (Mikkola & Aarseth 1996). The for-mer conserves the Hamiltonian and the angular momentum for a long-term evolution,while the latter significantly reduces the computational cost for a weakly perturbedbinary. In this work, the Hamilton equations of this algorithm are analyzed in detail.We mathematically and numerically show that it can correctly reproduce the secularevolution like the orbit averaged method and also well conserve the angular momen-tum. For a weakly perturbed binary, the method is possible to provide a few order ofmagnitude faster performance than the classical algorithm. A publicly available codewritten in the c++ language, sdar , is available on GitHub. It can be used either as astandalone tool or a library to be plugged in other N -body codes. The high precisionof the floating point to 62 digits is also supported. Key words: methods: numerical – software: simulations – Galaxy: globular clusters:general
Few-body dynamics embedded in N -body systems like starclusters is important for many research topics: the forma-tion of high-velocity stars via the ejections by interactingwith binaries or massive black holes (e.g. Poveda, Ruiz &Allen 1967; Leonard & Duncan 1988, 1990; Yu & Tremaine2003; Gualandris, Portegies Zwart & Sipior 2005; Gvara-madze, Gualandris & Portegies Zwart 2009; Fujii & Porte-gies Zwart 2011); the mergers of binaries that produce ex-otic 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 & Giersz2013) and gravitational wave sources (e.g. Portegies Zwart& McMillan 2000; O’Leary et al. 2006; Banerjee, Baum-gardt & Kroupa 2010; Antonini, et al. 2016; Askar, Szkud-larek, Gondek-Rosi´nska, Giersz & Bulik 2017; Di Carlo, etal. 2019); and the energy generation source that controls thedynamical evolution of star clusters (e.g. H´enon 1961; Heg-gie 1975; Hills 1975; Spitzer 1987; Binney & Tremaine 1987;Breen & Heggie 2013; Wang 2020). (cid:63) E-mail:[email protected]
However, few-body systems are also challenging to han-dle in the numerical simulations. Firstly, they can have amuch shorter dynamical timescale than that of the host en-vironment. For example, an open cluster with a few thou-sand stars has a typical dynamical (crossing) time of fewMyrs, but a close binary can have a period of few days. It isvery time consuming to accurately follow the lifetime of thecluster with a time resolution less than the binary period.Secondly, the numerical error introduced by the integratorcan accumulate after many binary orbits and cause an arti-ficial drift of the orbital parameters.To avoid these issues, the N -body codes for simulatingstar clusters, such as nbody6 (Aarseth 2003), artificially“freeze” the orbits of weakly perturbed binaries and hierar-chical systems, i.e., if the perturbation is below a thresholdand the hierarchical systems satisfy a stability criterion, theinternal motion of the systems are not evolved until the per-turbation becomes strong. Such trick can significantly re-duce the computational cost, but it ignores the long-termeffect of weak perturbation. In particular, the evolution ofangular momentum, and thus the eccentricity, can be com-pletely wrong. Besides, whether the instability of hierarchi- c (cid:13) a r X i v : . [ a s t r o - ph . E P ] F e b Long Wang et al. cal systems is properly treated depends on the quality of thestability criterion.There are two approximate solutions that can properlyhandle the few-body systems and also avoid the time con-suming calculation. The first method is using the orbit aver-aged Hamiltonian. In the case of a hierarchical triple, insteadof tracing the positions and velocities of individual compo-nents by integrating the classical Newtonian equation of mo-tion, the orbits of inner and outer binaries can be evolvedby using the Hamilton equation written in the Delaunay’selements. The short-period terms in the Hamiltonian can beeliminated (orbit averaged) by the Von Zeipel transforma-tion (e.g. Naoz, et al. 2013; Naoz 2016). For example, theorbit averaged method with a quadruple-level perturbationterm is provided in Naoz, et al. (2013). Because the timescaleof secular evolution is longer than the periods of binaries,such method can be very efficient to integrate the orbitalevolution of a stable hierarchical triple.Another method is called “slow-down” introduced byMikkola & Aarseth (1996). For a perturbed binary, by arti-ficially slowing down the orbital motion (scaling the time)while keeping the orbital parameters unchanged, the exter-nal perturbation is effectively enlarged. In such case, theeffect of perturbation on one binary orbit can represent theaverage effect of several orbits. When the perturbation isweak, this method can properly approximate the secular evo-lution. The benefit is rich: it is easy to be implemented in a N -body code; can be applied for a general few-body systemwith an arbitrary number of binaries and singles; is not onlyfor a kepler binary, but also for any type of periodic orbits.In Mikkola & Aarseth (1996), a constant slow-down factoris tested for a binary with a post-Newtonian type of pertur-bation. However, it is not well discussed and tested how wecan change the slow-down factor, even though it is crucialfor the actual use of the slow-down method. Especially, thisis important for dealing with the common few-body systemsin a star cluster, such as triples with high-eccentric or hy-perbolic orbits of perturbers. When the perturbation to abinary becomes smaller, we start from no slow-down to alarger slow-down factor. But how to change it properly isnot well understood. Mikkola & Aarseth (1996) discuss afactor-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 computingcost, but also avoid the large accumulate error due to muchless integration steps. But as the price of approximation,both methods lose the phase information of the orbits, thusthe mean motion resonance cannot be correctly reproduced.They also have different disadvantages. The orbit averagedmethod was used to study a limited type of stable hierar-chical systems, e.g. stable hierarchical triples and the Lunarsystem. For a more general case like the hyperbolic encoun-ters, systems with more than three bodies, or strongly per-turbed few-body systems, to derive the equation of motionwith high-order perturbation terms is unpractical. The slow-down method is much more flexible but it cannot avoid theissue of the singularity of Newtonian force, i.e., to integratethe motion of a high-eccentric binary, a large numerical er-ror appears at the peri-center unless the integration step issignificantly reduced.The solution for the issue of singularity is to apply theregularization algorithm. The key idea of regularization is to remove the singularity by a transformation of the equa-tion of motion. The Burdet-Heggie regularization (Burdet1967, 1968; Heggie 1973) for solving the Kepler problemuses a time-transformation, d t = r d τ , where τ is a ficti-tious time. By using the eccentric vector in the equationof motion together, the singularity term can be eliminated.Kustaanheimo & Stiefel (1965) introduce the KS regulariza-tion method, which transforms the equation of motion of theKepler orbit to a form of harmonic oscillator without singu-larity. Mikkola & Aarseth (1996) show how to combine theKS regularization together with the slow-down method tohave an efficient solution for integrating a few-body system.Mikkola & Tanikawa (1999) and Preto & Tremaine (1999)introduce the time-transformed explicit symplectic integra-tor (TSI) or “Algorithmic regularization” (AR), which canbe used for any potential that only depends on the coordi-nates of particles.Here we describe the TSI (AR) method in a bit moredetail. The symplectic integrator can conserve the Hamil-tonian and the angular momentum of a system. Thus it isvery suitable for simulating the long-term evolution of a sys-tem. However, it requires a constant integration step. In theclassical symplectic method, time step, d t , is also the inte-gration step, d s . This leads to a low efficiency in integratingan eccentric Kepler orbit. In order to be accurately enough,d t has to be fixed to the smallest value determined at thepericenter. The solution is to apply a time transformation,d t = g d s , which decouples d s and d t with a function g (e.gHairer 1997). Thus, d t can vary to avoid the issue of lowefficiency while d s is fixed to keep the symplectic property.This method can be described by the extended phase spaceHamiltonian, where t is treated as a coordinate and needto be integrated. The disadvantage is that the time trans-formation usually results in an inseparable Hamiltonian andonly the expensive implicit integrator can be used. Mikkola& Tanikawa (1999) and Preto & Tremaine (1999) find a so-lution by defining a specific type of g (see Section 3) thatthe Hamiltonian can be written in a separable style, in orderto use the explicit symplectic integrator. Especially, for anisolated binary system, if d t = d s/ | U | , where | U | is the ab-solute value of the binary potential energy (similar like theBurdet-Heggie time transformation), the integrator behav-iors dramatically well for the Kepler orbit, i.e., the numericaltrajectory follows the exact one with a phase error of time.In this work, we develop an efficient method for sim-ulating few-body systems by combining the slow-down andthe TSI (AR) schemes. In Section 2, we show the slow-downHamiltonian and the equation of motion for several typesof few-body systems, such as perturbed binaries, triples andgeneral hierarchical systems. We demonstrate that the slow-down method can correctly reproduce the secular evolutionand 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 dis-cussed in Section 4. The implementation and numerical testsare shown in Section 5 and 6. Finally, we discuss our resultsand draw conclusions in Section 7.To be convenient, hereafter r , p , v and m representcoordinates, conjugate momenta, velocities and masses ofparticles separately. We define the phase-space vector as w ≡ ( r , p ). The suffix of a variable using an index (e.g. i , j , 1, 2 ...) represents a specific particle while no suffix c (cid:13)000
However, few-body systems are also challenging to han-dle in the numerical simulations. Firstly, they can have amuch shorter dynamical timescale than that of the host en-vironment. For example, an open cluster with a few thou-sand stars has a typical dynamical (crossing) time of fewMyrs, but a close binary can have a period of few days. It isvery time consuming to accurately follow the lifetime of thecluster with a time resolution less than the binary period.Secondly, the numerical error introduced by the integratorcan accumulate after many binary orbits and cause an arti-ficial drift of the orbital parameters.To avoid these issues, the N -body codes for simulatingstar clusters, such as nbody6 (Aarseth 2003), artificially“freeze” the orbits of weakly perturbed binaries and hierar-chical systems, i.e., if the perturbation is below a thresholdand the hierarchical systems satisfy a stability criterion, theinternal motion of the systems are not evolved until the per-turbation becomes strong. Such trick can significantly re-duce the computational cost, but it ignores the long-termeffect of weak perturbation. In particular, the evolution ofangular momentum, and thus the eccentricity, can be com-pletely wrong. Besides, whether the instability of hierarchi- c (cid:13) a r X i v : . [ a s t r o - ph . E P ] F e b Long Wang et al. cal systems is properly treated depends on the quality of thestability criterion.There are two approximate solutions that can properlyhandle the few-body systems and also avoid the time con-suming calculation. The first method is using the orbit aver-aged Hamiltonian. In the case of a hierarchical triple, insteadof tracing the positions and velocities of individual compo-nents by integrating the classical Newtonian equation of mo-tion, the orbits of inner and outer binaries can be evolvedby using the Hamilton equation written in the Delaunay’selements. The short-period terms in the Hamiltonian can beeliminated (orbit averaged) by the Von Zeipel transforma-tion (e.g. Naoz, et al. 2013; Naoz 2016). For example, theorbit averaged method with a quadruple-level perturbationterm is provided in Naoz, et al. (2013). Because the timescaleof secular evolution is longer than the periods of binaries,such method can be very efficient to integrate the orbitalevolution of a stable hierarchical triple.Another method is called “slow-down” introduced byMikkola & Aarseth (1996). For a perturbed binary, by arti-ficially slowing down the orbital motion (scaling the time)while keeping the orbital parameters unchanged, the exter-nal perturbation is effectively enlarged. In such case, theeffect of perturbation on one binary orbit can represent theaverage effect of several orbits. When the perturbation isweak, this method can properly approximate the secular evo-lution. The benefit is rich: it is easy to be implemented in a N -body code; can be applied for a general few-body systemwith an arbitrary number of binaries and singles; is not onlyfor a kepler binary, but also for any type of periodic orbits.In Mikkola & Aarseth (1996), a constant slow-down factoris tested for a binary with a post-Newtonian type of pertur-bation. However, it is not well discussed and tested how wecan change the slow-down factor, even though it is crucialfor the actual use of the slow-down method. Especially, thisis important for dealing with the common few-body systemsin a star cluster, such as triples with high-eccentric or hy-perbolic orbits of perturbers. When the perturbation to abinary becomes smaller, we start from no slow-down to alarger slow-down factor. But how to change it properly isnot well understood. Mikkola & Aarseth (1996) discuss afactor-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 computingcost, but also avoid the large accumulate error due to muchless integration steps. But as the price of approximation,both methods lose the phase information of the orbits, thusthe mean motion resonance cannot be correctly reproduced.They also have different disadvantages. The orbit averagedmethod was used to study a limited type of stable hierar-chical systems, e.g. stable hierarchical triples and the Lunarsystem. For a more general case like the hyperbolic encoun-ters, systems with more than three bodies, or strongly per-turbed few-body systems, to derive the equation of motionwith high-order perturbation terms is unpractical. The slow-down method is much more flexible but it cannot avoid theissue of the singularity of Newtonian force, i.e., to integratethe motion of a high-eccentric binary, a large numerical er-ror appears at the peri-center unless the integration step issignificantly reduced.The solution for the issue of singularity is to apply theregularization algorithm. The key idea of regularization is to remove the singularity by a transformation of the equa-tion of motion. The Burdet-Heggie regularization (Burdet1967, 1968; Heggie 1973) for solving the Kepler problemuses a time-transformation, d t = r d τ , where τ is a ficti-tious time. By using the eccentric vector in the equationof motion together, the singularity term can be eliminated.Kustaanheimo & Stiefel (1965) introduce the KS regulariza-tion method, which transforms the equation of motion of theKepler orbit to a form of harmonic oscillator without singu-larity. Mikkola & Aarseth (1996) show how to combine theKS regularization together with the slow-down method tohave an efficient solution for integrating a few-body system.Mikkola & Tanikawa (1999) and Preto & Tremaine (1999)introduce the time-transformed explicit symplectic integra-tor (TSI) or “Algorithmic regularization” (AR), which canbe used for any potential that only depends on the coordi-nates of particles.Here we describe the TSI (AR) method in a bit moredetail. The symplectic integrator can conserve the Hamil-tonian and the angular momentum of a system. Thus it isvery suitable for simulating the long-term evolution of a sys-tem. However, it requires a constant integration step. In theclassical symplectic method, time step, d t , is also the inte-gration step, d s . This leads to a low efficiency in integratingan eccentric Kepler orbit. In order to be accurately enough,d t has to be fixed to the smallest value determined at thepericenter. The solution is to apply a time transformation,d t = g d s , which decouples d s and d t with a function g (e.gHairer 1997). Thus, d t can vary to avoid the issue of lowefficiency while d s is fixed to keep the symplectic property.This method can be described by the extended phase spaceHamiltonian, where t is treated as a coordinate and needto be integrated. The disadvantage is that the time trans-formation usually results in an inseparable Hamiltonian andonly the expensive implicit integrator can be used. Mikkola& Tanikawa (1999) and Preto & Tremaine (1999) find a so-lution by defining a specific type of g (see Section 3) thatthe Hamiltonian can be written in a separable style, in orderto use the explicit symplectic integrator. Especially, for anisolated binary system, if d t = d s/ | U | , where | U | is the ab-solute value of the binary potential energy (similar like theBurdet-Heggie time transformation), the integrator behav-iors dramatically well for the Kepler orbit, i.e., the numericaltrajectory follows the exact one with a phase error of time.In this work, we develop an efficient method for sim-ulating few-body systems by combining the slow-down andthe TSI (AR) schemes. In Section 2, we show the slow-downHamiltonian and the equation of motion for several typesof few-body systems, such as perturbed binaries, triples andgeneral hierarchical systems. We demonstrate that the slow-down method can correctly reproduce the secular evolutionand 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 dis-cussed in Section 4. The implementation and numerical testsare shown in Section 5 and 6. Finally, we discuss our resultsand draw conclusions in Section 7.To be convenient, hereafter r , p , v and m representcoordinates, conjugate momenta, velocities and masses ofparticles separately. We define the phase-space vector as w ≡ ( r , p ). The suffix of a variable using an index (e.g. i , j , 1, 2 ...) represents a specific particle while no suffix c (cid:13)000 , 1–14 DAR integrator represents all particles in a N -body system. For exmample, w ≡ ( w , w ,..., w N ). The slow-down method can be described by a modifiedHamiltonian (Mikkola & Aarseth 1996). For a perturbedbinary, the slow-down Hamiltonian is H sd = 1 κ H b + ( H − H b ) , (1)where H is the original Hamiltonian of the system, H b is theHamiltonian of the binary and κ is the slow-down factor. Theequation of motion is d w d t = { w , H sd } , (2)where {} is the Possion bracket. When κ = 1, H sd ≡ H .The important point is that although the Hamiltonian ismodified, w keeps the original definition. Thus, the slow-down affects the time evolution of orbital phase, but thepositions and velocities of particles are not scaled. We first consider a simple example of one binary with anexternal potential U ( r ) and a fixed κ , H sd = 1 κ (cid:20) m v + 12 m v − Gm m | r − r | (cid:21) + U ( r ) . (3)where the suffixes “1” and “2” denote the two componentsof the binary, and p i is replaced by m i v i .Applying Eq. 2, the evolution of r and v can be de-scribed byd v i d t = − κ Gm − i ( r i − r − i ) | r i − r − i | − ∂U ( r ) ∂ r i ( i = 1 or 2)d r i d t = 1 κ v i . (4)Thus, in the reference frame of the perturber, the motion ofbinary is slowed down by a factor of κ , while in the view ofthe binary, the perturbation is enhanced by κ times. We canunderstand this better by considering the cumulative effectof perturbation to v i on one slow-down orbit.The effective period of the slow-down binary has P sd = κP , where P is the original period. By integrating d v i / d t ofEq. 4, the change of v i after one P sd is δ v i = (cid:90) P sd d v i d t . (5)In the case of no perturbation, v i ( t = P sd ) = v i ( t = 0).Thus (cid:90) P sd d v i d t = 0 . (6)When the effect of U ( r ) is weak, δ v i ≈ (cid:90) P sd − ∂U ( r ) ∂ r i . (7)Then δ v i is only determined by the integrated effect of per-turbation. Here r i only passes one cycle, but the integration time is κP . Thus, the effect the perturbation is amplified by κ times on one orbit of binary.On the other hand, when κ is an integer, we can mea-sure δ v i in the original case (no slow-down; denoted as δ v (cid:48) i )after the same physical time by a summation of κ pieces ofintegral, where each piece represents one orbit: δ v (cid:48) i = κ (cid:88) j =1 (cid:90) jP ( j − P − ∂U ( r (cid:48) ) ∂ r (cid:48) i = κ (cid:88) j =1 V j , (8)If the perturbation does not change significantly, V j ≈ V ( j = 1 , κ ) . (9)This means the effect of perturbation can be representedby the integration of one orbit ( V ) multiplied by κ times.Thus, it is equivalent to the treatment of the slow-downmethod and δ v i ≈ δ v (cid:48) i . In other words, the slow-downmethod approximates the perturbation in an orbit averagedway. However, when the external potential changes signifi-cantly within the time interval of κP , Eq. 9 becomes invalid.Thus, there is a threshold of κ . We discuss the criterion inSection 2.4. For a hierarchical triple, H sd = 1 κ (cid:20) m ( v − v cm ) + 12 m ( v − v cm ) − Gm m | r − r | (cid:21) + 12 ( m + m ) v + 12 m v − Gm m | r − r | − Gm m | r − r | , (10)where the suffixes “1” and “2” denote the two componentsof the inner binary and “3” represents the third body, and v cm is the center-of-the-mass velocity of the inner binary.Because the slow-down affects only the internal motion ofthe binary in its rest frame, v cm is subtracted from the slow-down term.Applying Eq. 2 (with a fixed κ to Eq. 10), the evolutionof r and v can be described by • binary components ( i = 1 or 2):d v i d t = − κ Gm − i ( r i − r − i ) | r i − r − i | − Gm ( r i − r ) | r i − r | d r i d t = 1 κ ( v i − v cm ) + v cm ; (11) • third body ( i = 3):d v i d t = − (cid:88) j =1 Gm j ( r i − r j ) | r i − r j | d r i d t = v i . (12)The form of the third body is identical to the original one.We demonstrate how the perturbation force is calcu-lated in the slow-down and the original integration methodsin Fig. 1. Along the binary orbit (one of the component), 8sample points with an equal angular separation are chosenfor the demonstration. When the binary finish two cycles c (cid:13) , 1–14 Long Wang et al.
123 4 5 6 78
Figure 1.
The illustration to show in a hierarchical triple howthe force from the perturber to the binary component is calcu-lated in the slow-down and the original methods. The left circlerepresents the orbit of one binary component and the right curveindicates the track of the perturber. The points alone the twotracks indicate the corresponding positions (same color or num-ber) for calculating the perturbation force. The number representsthe original method where binary pass two cycles while the colorrepresents the slow-dowm method with κ = 2 (one cycle). Forexample, the two solid lines and one dashed line link the cor-responding positions (dark blue point with number of 7) in theoriginal and the slow-down cases separately. in the original case, each sample point has twice force cal-culations from the perturber at the corresponding positionsshown with the same number along its track. In the slow-down case with κ = 2, the binary finish one orbit when theperturber passes the same track. Thus, the correspondingposition of the perturber (with same colors) for each samplepoint is different from that of the original case. If the samplepoints indicate the step sizes in an integrator, the slow-downmethod reduces the step sizes by κ times.The two types of lines that connect the dark blue point7 represent one example of the perturbation calculation inthe two methods separately. If the solid and dashed lineshave very close length and direction, the slow-down methodgive a good approximation of the perturbation force. This isthe case when the perturber is far away ( | r | (cid:29) | r − r | ),i.e., the perturbation is weak. For a hierarchical triple, we can also describe the slow-downmethod in terms of the orbit averaged method. By using theLegendre expansion for the perturbation term (e.g. Harring-ton 1968; Naoz, et al. 2013; Naoz 2016), H sd can be writtenas the combination of three components: H sd = 1 κ H b , in + H b , out + H pert H b , in = − Gm m a in H b , out = − G ( m + m ) m a out H pert = − G | r out | ∞ (cid:88) n =2 M n (cid:18) | r in || r out | (cid:19) n P n (cos Ψ) , (13) where P n is the Legendre polynomials, the suffixes “in” and“out” indicate the inner and outer binaries, a represents thesemi-major axis, r in = r − r , r out = r − r cm , r cm is thecenter-of-the-mass position of the inner binary and Ψ is theangle between r in and r out . The mass coefficient sequence, M n = m m m m n − − ( − m ) n − ( m + m ) n . (14)Using the Delaunay’s elements, the two binary terms ofHamiltonian have the form H b , in = − G m m m + m ) L H b , out = − G ( m + m ) m m + m + m ) L (15)where the conjugate momenta, L in = m m m + m (cid:112) G ( m + m ) a in L out = ( m + m ) m m + m + m (cid:112) G ( m + m + m ) a out . (16)Applying the equation of motion, the time derivations oftheir corresponding coordinates, (cid:96) in and (cid:96) out , ared (cid:96) in d t = ∂H sd ∂ L in = 1 κ G m m ( m + m ) L + ∂H pert ∂ L in d (cid:96) out d t = ∂H sd ∂ L out = G ( m + m ) m ( m + m + m ) L + ∂H pert ∂ L out . (17)The effect of slow-down is reflected on the first termof d (cid:96) in / d t , which represents the mean motion of the innerbinary. Without the perturbation term, Eq. 17 indicates thatthe mean motion of the inner binary is κ times slower ofthat in the original case, as discussed in Section 2.1. Theorbit averaged method uses a canonical transformation (VonZeipel transformation) to remove the short-period terms ( (cid:96) in and (cid:96) out ) in H pert (e.g. Naoz, et al. 2013): H pert = (cid:90) π (cid:90) π H pert d (cid:96) in d (cid:96) out . (18)Since (cid:96) in and (cid:96) out do not appear in H b , in and H b , out and H pert is independent of (the fixed) κ , the transformationresults in the same form of H pert in the slow-down and orig-inal cases. Thus, the slow-down method provides the samesecular evolution as the orbit averaged way. N -body system with multiple slow-downbinaries For a general few-body system containing N b binaries, eachbinary can have an individual slow-down factor, κ i . Theslow-down Hamiltonian has the general form: H sd = N b (cid:88) i κ i H b ,i + (cid:32) H − (cid:88) i H b ,i (cid:33) H b ,i = 12 m ,i ( v ,i − v cm ,i ) + 12 m ,i ( v ,i − v cm ,i ) − Gm ,i m ,i | r ,i − r ,i | v cm ,i = m ,i v ,i + m ,i v ,i m ,i + m ,i , (19) c (cid:13)000
The illustration to show in a hierarchical triple howthe force from the perturber to the binary component is calcu-lated in the slow-down and the original methods. The left circlerepresents the orbit of one binary component and the right curveindicates the track of the perturber. The points alone the twotracks indicate the corresponding positions (same color or num-ber) for calculating the perturbation force. The number representsthe original method where binary pass two cycles while the colorrepresents the slow-dowm method with κ = 2 (one cycle). Forexample, the two solid lines and one dashed line link the cor-responding positions (dark blue point with number of 7) in theoriginal and the slow-down cases separately. in the original case, each sample point has twice force cal-culations from the perturber at the corresponding positionsshown with the same number along its track. In the slow-down case with κ = 2, the binary finish one orbit when theperturber passes the same track. Thus, the correspondingposition of the perturber (with same colors) for each samplepoint is different from that of the original case. If the samplepoints indicate the step sizes in an integrator, the slow-downmethod reduces the step sizes by κ times.The two types of lines that connect the dark blue point7 represent one example of the perturbation calculation inthe two methods separately. If the solid and dashed lineshave very close length and direction, the slow-down methodgive a good approximation of the perturbation force. This isthe case when the perturber is far away ( | r | (cid:29) | r − r | ),i.e., the perturbation is weak. For a hierarchical triple, we can also describe the slow-downmethod in terms of the orbit averaged method. By using theLegendre expansion for the perturbation term (e.g. Harring-ton 1968; Naoz, et al. 2013; Naoz 2016), H sd can be writtenas the combination of three components: H sd = 1 κ H b , in + H b , out + H pert H b , in = − Gm m a in H b , out = − G ( m + m ) m a out H pert = − G | r out | ∞ (cid:88) n =2 M n (cid:18) | r in || r out | (cid:19) n P n (cos Ψ) , (13) where P n is the Legendre polynomials, the suffixes “in” and“out” indicate the inner and outer binaries, a represents thesemi-major axis, r in = r − r , r out = r − r cm , r cm is thecenter-of-the-mass position of the inner binary and Ψ is theangle between r in and r out . The mass coefficient sequence, M n = m m m m n − − ( − m ) n − ( m + m ) n . (14)Using the Delaunay’s elements, the two binary terms ofHamiltonian have the form H b , in = − G m m m + m ) L H b , out = − G ( m + m ) m m + m + m ) L (15)where the conjugate momenta, L in = m m m + m (cid:112) G ( m + m ) a in L out = ( m + m ) m m + m + m (cid:112) G ( m + m + m ) a out . (16)Applying the equation of motion, the time derivations oftheir corresponding coordinates, (cid:96) in and (cid:96) out , ared (cid:96) in d t = ∂H sd ∂ L in = 1 κ G m m ( m + m ) L + ∂H pert ∂ L in d (cid:96) out d t = ∂H sd ∂ L out = G ( m + m ) m ( m + m + m ) L + ∂H pert ∂ L out . (17)The effect of slow-down is reflected on the first termof d (cid:96) in / d t , which represents the mean motion of the innerbinary. Without the perturbation term, Eq. 17 indicates thatthe mean motion of the inner binary is κ times slower ofthat in the original case, as discussed in Section 2.1. Theorbit averaged method uses a canonical transformation (VonZeipel transformation) to remove the short-period terms ( (cid:96) in and (cid:96) out ) in H pert (e.g. Naoz, et al. 2013): H pert = (cid:90) π (cid:90) π H pert d (cid:96) in d (cid:96) out . (18)Since (cid:96) in and (cid:96) out do not appear in H b , in and H b , out and H pert is independent of (the fixed) κ , the transformationresults in the same form of H pert in the slow-down and orig-inal cases. Thus, the slow-down method provides the samesecular evolution as the orbit averaged way. N -body system with multiple slow-downbinaries For a general few-body system containing N b binaries, eachbinary can have an individual slow-down factor, κ i . Theslow-down Hamiltonian has the general form: H sd = N b (cid:88) i κ i H b ,i + (cid:32) H − (cid:88) i H b ,i (cid:33) H b ,i = 12 m ,i ( v ,i − v cm ,i ) + 12 m ,i ( v ,i − v cm ,i ) − Gm ,i m ,i | r ,i − r ,i | v cm ,i = m ,i v ,i + m ,i v ,i m ,i + m ,i , (19) c (cid:13)000 , 1–14 DAR integrator where the suffixes “1 , i ” and “2 , i ” represent the two com-ponents in the i th binary. To be convenient, we use the first2 N b indices to indicate binary components and others to in-dicate singles, i.e., “1,i” and “2,i” are equivalent to 2 i − i separately.Applying Eq. 2 with all κ i being fixed, the equation ofmotion is • binary components ( i = 1 , , ..., N b ; k = 1 or 2):d v k,i d t = − κ i Gm − k,i ( r k,i − r − k,i ) | r k,i − r − k,i | − N (cid:88) j =1; j (cid:54) =2 i +1 − k Gm j ( r k,i − r j ) | r k,i − r j | d r k,i d t = 1 κ i ( v k,i − v cm ,i ) + v cm ,i ; (20) • single body ( i = 2 N b + 1 , N b + 2 , ..., N ):d v i d t = − N (cid:88) j =1; j (cid:54) = i Gm j ( r i − r j ) | r i − r j | d r i d t = v i . (21)This general form suggests that κ i only explicitly affects theinternal motion of the binary i . The force calculations ofsingles and other binaries do not explicitly depend on κ i . κ in the integration In previous sections, we assume that κ is constant. In manyfew-body systems, the perturbation to a binary can vary sig-nificantly, especially when a nearby perturber has a highlyeccentric or hyperbolic orbit. Thus, varying κ during theintegration is necessary to balance the performance and ac-curacy. Mikkola & Aarseth (1996) suggests to calculate κ based on the ratio between the tidal force from perturbersand the internal force of the binary assuming the separationis 2 a in . For a system with one binary and N p perturbers, wedefine a modified version by including the eccentricity ( e ), κ ( w ) = k ref m m ( m + m ) [ a in (1 + e in )] N p (cid:88) i | r i − r cm | m i , (22)where k ref is a constant coefficient.Since this κ depends on w , one might think that the ad-ditional terms of ( H b , in ∂κ ( w ) /∂ r , H b , in ∂κ ( w ) /∂ v ) shouldappear in the equation of motion. However, these terms ac-tually should not be included to ensure a correct secularevolution. In the case of a triple system, by applying Eq. 22in the Hamiltonian form of Eq. 13,1 κ H b , in = − k ref G ( m + m ) m | r out | (cid:20) a in (1 + e in ) | r out | (cid:21) . (23) H b , in becomes to depend on (cid:96) out (from a in / | r out | ) in thesame order as that of the perturbation term in H pert with n = 2 . The additional terms including H b , in then affectsthe orbital average of (cid:96) out (Eq. 18), which may not result inthe same secular evolution as the orbit averaged method.On the other hand, when Eq. 9 is satisfied, κ ( w ) isexpected to evolve slowly in one step. As the first order ap-proximation, κ can be treated as a temporary constant and isonly updated at the end of one integration step. Essentially, κ is only a step function of time. Then the unwarranted ef-fect can be avoided. We can describe this Jumping- κ methodfor a system with N b binaries as:(i) advance d t with fixed κ i measured at t .(ii) update κ i at the new time t (cid:48) = t + d t ;(iii) apply an instant correction of H sd by δH sd ( t (cid:48) ) = N b (cid:88) i (cid:18) κ i ( t (cid:48) ) − κ i ( t ) (cid:19) H b ,i ( t (cid:48) ) . (24)With the step (iii), even H sd is not conserved, the cumulativenumerical error of H sd can still be properly tracked as (cid:15) ( H sd , k d t ) = H sd ( k d t ) − k (cid:88) i =1 δH sd ( i d t ) − H sd (0) , (25)where k is the total step count.In this method, although the unwarranted effect on thesecular evolution is not included in one integration step, theinstant change of H sd may still introduce an additional ef-fect. Thus, it is necessary to limit the change rate of κ forsafety. We estimate how κ ( w ) changes depending on the or-bital parameters. With Eq. 22, the differential of 1 /κ ,d (cid:18) κ (cid:19) = k ref N p (cid:88) i m + m ) m i m m a | r i − r cm | (cid:18) d aa − d | r i − r cm || r i − r cm | (cid:19) . (26)Typically, in the weak perturbation condition, a does notchange significantly but | r i − r cm | may vary a lot if theperturber has an eccentric or hyperbolic orbit. In the caseof one perturber, Eq. 26 can be simplified asd (cid:18) κ (cid:19) ≈ − k ref (cid:18) κ (cid:19) d | r − r cm || r − r cm | . (27)Thus, the time derivation of 1 /κ depends on | v − v cm | / | r − r cm | . The evolution timescale is the order of the crossingtime of the perturber.The slow change of κ requires that the strength of per-turbation does not vary significantly within at least one P sd .Based on Eq. 27, we can introduce a timescale criterion tolimit the maximum value of κ for safety: κ max = c |(cid:104) r i − r (cid:105)| P sd |(cid:104) v i − v cm (cid:105)| , (28)where c is a constant coefficient; (cid:104) r i − r (cid:105) and (cid:104) v i − v cm (cid:105) represent the mass weighted average of velocities and posi-tions of all perturbers referring to the center-of-the-mass ofthe binary separately. If a binary is the perturber, the massweights help to remove the velocity fluctuations caused bythe internal motion of the two components of the perturber.For the Jumping- κ method, Eq. 22 and 28 together providethe criterion to estimate κ . In Section 6, we show that thenumerical tests by using the Jumping- κ method indeed pro-vide a correct secular evolution. It is important to know, with the irregular form of H sd ,what are conserved quantities during the evolution that can c (cid:13) , 1–14 Long Wang et al. be used to check the quality of the integration. When κ is a constant, H sd does not explicitly depend on time. Thusthe “slow-down energy”, H sd , is a conserved quantity, whichmeans that the physical energy, H , changes during the evo-lution. The variation of H can be calculated by δH = H − H sd = N b (cid:88) i (cid:18) − κ i (cid:19) H b ,i . (29)Thus, H oscillates in a timescale of the secular evolution ofthe binaries.When κ is not fixed, whether H sd is conserved de-pends on how κ is evaluated. If κ follows Eq. 22 exactlyand the equation of motion contains the additional termsof ( H b , in ∂κ ( w ) /∂ r , H b , in ∂κ ( w ) /∂ v ), κ would not explic-itly depend on time, thus H sd is conserved. However, thisalso indicates that a significant variation of κ results in alarge change of binary energy. Such behaviour is introducedby the unwarranted effect discussed in Section 2.4. Instead,the Jumping- κ method that exclude the additional term canavoid such problem, although H sd explicitly depends on timeand is not conserved. Since the physical energy is anywaynot conserved, it is more important to ensure the correctbehaviour of secular evolution rather than to conserve H sd .On the other hand, we can still follow the integration errorby using Eq. 25 in the absence of the energy conservation. The angular momentum defined by L = (cid:88) i m i r i × v i (30)is a conserved quantity in the original case since { L , H } = , (31)where is a zero vector. Interestingly, it can be proved that L is also conserved with H sd (Eq. 19; see Appendix), i.e. { L , H sd } = . (32)This conservation is not only valid for an invariant κ . Withthe help of a symbolic calculator , it can be shown that for κ defined in Eq. 22, Eq. 32 is also satisfied. On the other hand,if κ explicitly depends on time, L is still conserved becausethe differential in Eq. 32 is independent of time. Thus, L is generally conserved in the slow-down method and can beused to validate the quality of the integration. The TSI (AR) method is accurate for solving the long-termevolution of few-body systems. It can be described by theextended phase-space Hamiltonian (e.g. Hairer 1997; Preto& Tremaine 1999):Γ( W ) = g ( W )[ H ( w , t ) − H ( w (0) , , (33)where H ( w , t ) is the standard Hamiltonian and g ( W ) istime-transformation function. The extended phase-space We use python.sympy to prove Eq. 32. * Figure 2.
The schedule of the LogH method with the equationof motion shown in Eq. 38 and 39. vector, W = ( R , P ), is defined by including an additionalpair of the coordinate, t , and the conjugate momentum, p t = − H ( w (0) ,
0) (negative initial energy). w (0) is the ini-tial value of w .With the new differential variable s , g ( W ) = d t d s . (34)A special type of g ( W ) introduced by Preto & Tremaine(1999), g ( W ) = f ( T ( P )) − f ( − U ( R )) T ( P ) + U ( R ) , (35)leads to a separable Γ:Γ( W ) = f ( T ( P )) − f ( − U ( R )) , (36)where the kinetic energy, T ( P ) ≡ T ( p ) + p t , and the poten-tial energy, U ( R ) ≡ U ( r , t ).Putting Γ( W ) in the equation of motion,d W d s = { W , Γ( W ) } , (37)the derivatives of W with respect to s ared r i d s = f (cid:48) ( T ( p ) + p t ) ∂T ( p ) ∂ p i d t d s = f (cid:48) ( T ( p ) + p t )d p i d s = f (cid:48) ( − U ( r , t )) ∂U ( r , t ) ∂ r i d p t d s = f (cid:48) ( − U ( r , t )) ∂U ( r , t ) ∂t . (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 canensure that the numerical trajectory of a Kepler orbit fol-lows the exact one, w ( t ), with only a phase error of time.Hereafter we use “LogH” to represent this algorithm. If thetime error is not important, the LogH method is very effi-cient to integrate a weakly perturbed Kepler orbit. Fig. 2shows the schedule of one DKD loop.When H does not explicitly depends on t , H ( w , t ) = c (cid:13)000
0) (negative initial energy). w (0) is the ini-tial value of w .With the new differential variable s , g ( W ) = d t d s . (34)A special type of g ( W ) introduced by Preto & Tremaine(1999), g ( W ) = f ( T ( P )) − f ( − U ( R )) T ( P ) + U ( R ) , (35)leads to a separable Γ:Γ( W ) = f ( T ( P )) − f ( − U ( R )) , (36)where the kinetic energy, T ( P ) ≡ T ( p ) + p t , and the poten-tial energy, U ( R ) ≡ U ( r , t ).Putting Γ( W ) in the equation of motion,d W d s = { W , Γ( W ) } , (37)the derivatives of W with respect to s ared r i d s = f (cid:48) ( T ( p ) + p t ) ∂T ( p ) ∂ p i d t d s = f (cid:48) ( T ( p ) + p t )d p i d s = f (cid:48) ( − U ( r , t )) ∂U ( r , t ) ∂ r i d p t d s = f (cid:48) ( − U ( r , t )) ∂U ( r , t ) ∂t . (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 canensure that the numerical trajectory of a Kepler orbit fol-lows the exact one, w ( t ), with only a phase error of time.Hereafter we use “LogH” to represent this algorithm. If thetime error is not important, the LogH method is very effi-cient to integrate a weakly perturbed Kepler orbit. Fig. 2shows the schedule of one DKD loop.When H does not explicitly depends on t , H ( w , t ) = c (cid:13)000 , 1–14 DAR integrator Figure 3.
The schedule of the Mikkola (2012) scheme. H ( w (0) , T ( P ) + U ( R ) = 0. Thus f (cid:48) ( T ( P )) = f (cid:48) ( − U ( R )). Mikkola & Aarseth (2002) found that insteadof calculating T ( P ), one can also define a variable, u = (cid:90) ∂U ( R ) ∂ R · d R d t (40)Then f (cid:48) ( u ) = f (cid:48) ( T ( P )). The integration schedule of thismethod for one DKD step is shown in Fig. 3. Via the Taylor expansion, Eq. 35 can be approximated as(Preto & Tremaine 1999) g ( W ) ≈ f (cid:48) ( − U ( R )) (41)By applying Eq. 39 for a Kepler orbit system, the relationbetween t and s isd t ≈ d s − U ( R ) = | r − r | d sGm m . (42)This is equivalent to the case in Burdet-Heggie regulariza-tion method (Burdet 1967, 1968; Heggie 1973). It can beshown that this time transformation can remove the singu-larity by including the energy and the eccentric vector in theequation of motion. This suggests that the LogH method canwell handle the high-eccentric orbit.On the other hand, we can show that s has a geometricmeaning. For a Kepler orbit. | r − r | = 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 = P π ( E − e sin E ) . (44)The differential of this equation results ind t = P π (1 − e cos E )d E. (45)By using Eq. 43, the time derivative of E isd E d t = a | r − r | πP . (46)This suggests that s represents a scaled E :d s = m m (cid:114) Gam + m d E = L d E, (47) Thus, the fixed d s in the symplectic integrator indicates aconstant step of d E . Since t is unknown before the integra-tion finish, it is difficult to determine d s to stop the inte-gration at a given t . However, the geometric meaning of s provides a way to stop at a certain E . This is very usefulfor determining the number of steps per orbit in order tocontrol the numerical accuracy of the integration. Although the LogH method ensures that the numerical tra-jectory follows the exact Kepler orbit, the numerical H ( w )is not perfectly conserved. Especially, when the orbit reachesthe pericenter, the error of H ( w ) is large due to the trunca-tion of floating points. Actually, based on the definition, theexact conserved “energy” is Γ( W ). With Eq. 33, 39 and 41,Γ( W ) ≈ | r − r | Gm m [ H ( w , t ) − H ( w (0) , . (48)Since | r − r | reaches the minimum at the pericenter, thelarge error of H ( w ) is cancelled out in Γ( W ). In Section 6,we show the numerical test of such behaviour. This featureindicates that we should avoid the calculation the orbitalparameters or stop the integration at the pericenter in orderto avoid a large numerical error of the orbital elements. The slow-down method helps to solve the large timescale is-sue of a hierarchical system, while the TSI (AR) method canconserve Γ and 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,Γ( W ) can take any separable Hamiltonian, thus it is straightforward to implement H sd ( w , t ) into Γ( W ). With Eq. 19 and36, we can obtain the time transformed slow-down Hamil-tonian for a N -body system with N b binaries:Γ sd ( W ) = f ( T sd ( P )) − f ( − U sd ( R )) T sd ( P ) = N b (cid:88) i =1 κ i T b ,i ( v ,i , v ,i ) + N b (cid:88) i =1 T b , cm ,i ( v cm ,i )+ N (cid:88) i =2 N b +1 T i ( v i ) + p t U sd ( P ) = N b (cid:88) i =1 κ i U b ,i ( r ,i , r ,i ) + (cid:34) U ( r ) − N b (cid:88) i = i U b ,i ( r ,i , r ,i ) (cid:35) , (49)where T b ,i ( v ,i , v ,i ) = 12 m ,i ( v ,i − v cm ,i ) + 12 m ,i ( v ,i − v cm ,i ) T b , cm ,i ( p cm ,i ) = 12 ( m ,i + m ,i ) v ,i U b ,i ( r ,i , r ,i ) = − Gm ,i m ,i | r ,i − r ,i | U ( r ) = − N (cid:88) i =1 N (cid:88) j = i +1 Gm i m j | r i − r j | . (50) c (cid:13) , 1–14 Long Wang et al.
The equation of motion of Γ sd ( W ) is the same as Eq. 38except that T ( P ) and U ( P ) are replaced by T sd ( P ) and U sd ( P ) separately.Notice that here the conserved “energy” is Γ sd (withfixed κ ). Thus when the Jumping- κ method (see Section 2.4)is applied, the correction term in the third step should bemodified as δ Γ sd ( t (cid:48) ) = f (cid:2) T sd [ P ( t (cid:48) ) , κ ( t (cid:48) )] (cid:3) − f (cid:2) U sd [ R ( t (cid:48) ) , κ ( t (cid:48) )] (cid:3) − (cid:2) f (cid:2) T sd [ P ( t (cid:48) ) , κ ( t )] (cid:3) − f (cid:2) U sd [ R ( t (cid:48) ) , κ ( t )] (cid:3)(cid:3) (51)The corresponding cumulative numerical error, ε (Γ sd ), canbe obtained by replacing H sd to Γ sd in Eq. 25. Based on the SDAR method, we develop a software library, sdar , written in the C++ language with the object orientedprogramming . This library contains three modules: ar , hermite and kepler . ar is the implementation of the SDAR method shownin Eq. 50 with the Jumping- κ algorithm. The high-orderexplicit 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 κ i ,while the system as a whole can also apply the slow-downmethod if it has a periodical evolution. Thus, hermite al-lows a nested two-level slow-down treatment. This can ac-celerate the integration of a weakly perturbed stable hierar-chical system. kepler is a tool to construct a Hierarchical (Kepler)binary tree for a group of particles. It is used to identify theinner binaries and calculate κ .The quadruple-double precision library, qd (Hida, Li& Bailey 2001), is included in the code. Thus the user canswitch on the high-precision support (up to 62-digit pre-cision). Notice that when qd is used, the performance isreduced and also it is not thread-safe. In this section, we check whether the slow-down methodcan correctly reproduce the secular evolution of few-bodysystems. The comparison between the slow-down and orig-inal algorithms is done by using the ar module with the6 th -order symplectic method (Solution A in Table 1 fromYoshida 1990). We use the Delaunay’s elements to describethe geometric information of the binary orbits. By selectingthe Cartesian coordinate system ( x - y - z ), the three anglesare (see Fig. 4): • θ : inclination; the angle from z -axis to L . URL: https://github.com/lwang-astro/SDAR x yz x' L θ φ ψ e Figure 4.
The Delaunay’s elements (three Euler angles) of abinary orbit in the Cartesian coordinate system ( x - y - z ). θ : incli-nation, φ : longitude of the ascending node, and ψ : argument ofperiapsis. • φ : longitude of the ascending node (often denoted asΩ); the angle from x -axis to the intersection of the orbitalplane and the x - y plane ( x (cid:48) -axis). • ψ : argument of periapsis (often denoted as ω ); the anglefrom eccentric vector to the x (cid:48) -axis. The Kozai-Lidov effect is one of the most important featureappearing in a family of hierarchical triple systems (Kozai1962; Lidov 1962). For example, when the outer binary hasa circular orbit and the secondary of the inner binary hasa negligible mass compared to the other two bodies, the z (cid:48) component of the orbit averaged angular momentum of theinner binary is a conserved quantity: L z (cid:48) , ∝ (cid:113) (1 − e ) cos( θ in ) = const , (52)Where the z axis is defined in the invariable plane. When θ in is in the range of cos − ( ± (cid:112) / e in and θ in can happen. In the astrophysical environment,such oscillation can result in a high eccentricity that maytrigger the merger of the inner binary. Thus, it is impor-tant to validate that the slow-down method can correctlyreproduce the Kozai-Lidov effect.We perform a simulation of a hierarchical triple (B-S).The formula from Antognini (2015) is used to estimate theKozai-Lidov oscillation timescale: t KL ≈ π (cid:18) m m (cid:19) P ut P i n (1 − e ) / . (53)For the purpose of test, parameters which allow a large κ anda short t KL is preferred. Thus, smaller m /m and P o ut /P i n or higher eccentricity of outer binary are better. However,the former also reduce the maximum κ (Eq. 22), thus thelatter is good for both requirements. c (cid:13) , 1–14 DAR integrator a ×10 +10 in SD ORG0.9991.0001.001 out e
101 ×10 +9.9×10 Figure 5.
The evolution of orbital parameters of the inner andouter binaries for the B-S system. The red color represent theresult using the slow-down method (SD) and the black color rep-resents the case of no slow-down (ORG). For each panel, we applythe scientific notation in the plotting style of y -axis: the actualvalues of y -axis are calculated by y tick × scale + y offset , where y tick is the value shown along the the y-axis, scale is the firstvalue shown above the y-axis ( scale = 1 in default) and y offset isthe second value following the symbol “+” ( y offset = 0 in default). One suitable initial condition is shown in Table 1, where t KL ≈
45. For convenience, we use the scale-free unit (gravi-tational constant is one). This is also applied for all numer-ical tests discussed below.Both the inner and outer binaries are initially near theapo-centre positions ( E ≈ . k ref = 1 . × − as suggested by Mikkola & Aarseth (1996) in Eq. 22, κ max ≈
52 when the outer body is at the apo-centers. Asthe outer eccentricity, e out , is large, κ decreases to belowone based on Eq. 22 at the pericenter position. This shouldbe avoided, thus we let κ (cid:62) .
0. To obtain a high numericalaccuracy, we use 256 steps per orbit of the inner binary.This results in ∆ s ≈ . × − . To reach about 4 t KL , thetotal number of integration step without slow-down is about2 . × . In the slow-down case, it is about 3 . × stepswhich is 6 times smaller.Fig. 5 shows the simulation result of the B-S systemby using the LogH method with and without slow-down(hereafter named as “SD” and “ORG” methods) . TheKozai-Lidov oscillation appears in the evolution of the or-bital parameters. The comparison between red (SD) andblack (ORG) curves clearly indicates that the SD methodcan well reproduce the Kozai-Lidov effect with the correcttimescale. All orbital parameters, including a , e and three The scientific notation in the plotting style of y -axis (describedin Fig. 5) is defined by the python module matplotlib.ticker (see https://matplotlib.org/ for reference). In all figures, we followthis style in order to save space. a ×10 +9.999999×10 e ×10 +9×10 Figure 6.
The evolution of a and e for the inner binary of the B-Ssystem with a high-time-resolution by using the ORG method. Euler angles of both inner and outer binaries are consistentwith each others. The oscillation timescale is not exacted aspredicted by Eq. 53. Since the secondary of the inner binaryhas a small but none zero mass, the B-S system does notexactly satisfy the test-particle condition that is assumed inEq. 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. 6shows the high-time-resolution evolution of a in and e in with t in 0 ∼ × − by using the ORG method. When theinner binary passes the pericenter, sharp peaks appear. Theamplitude of peaks depends on the orbital phase of bothinner and outer binaries. Since the time interval of one peakis very short, most of the peaks are not sampled in the low-time-resolution plot and they appear by chance. These peaksalso exist in other orbital parameters but cannot be seen inthe plots due to a large amplitude of the secular variation.As the peaks do not affect the secular evolution, they canbe safely ignored. κ oscillates when the third body cycles between the apo-center and the pericenter, as shown in the upper panel ofFig. 7. The evolution of H sd follows the pattern of κ asdescribed by Eq. 29. Although the scatter of H sd is large, thebehaviour of the secular evolution is correct. This indicatesthat whether H sd is conserved does not matter.However, we still can trace the numerical error by us-ing ε (Γ sd ). In Fig. 8 the evolution of the cumulative nu-merical errors of different definitions of Hamiltonian, ε ( H ), ε ( H sd ) and ε (Γ sd ), are compared. ε ( H ) is the error of origi-nal Hamiltonian, defined as H ( t ) − H (0). After e in pass themaximum value for the first time, large oscillations appearson both ε ( H ) and ε ( H sd ). ε ( H ) of the SD method has amuch larger scatter compared to the ORG method, while inthe case of ε ( H sd ), both have a similar amplitude. This isconsistent with the expectation as discussed in Section 2.5.1.On the other hand, the scatter of ε (Γ) or ε (Γ sd ) is two or-der of magnitude smaller than ε ( H sd ) without oscillations,which indicates that the conserved quantities are Γ (ORG) c (cid:13) , 1–14 Long Wang et al.
Table 1.
The initial condition of the hierarchical triple (B-S) system for test. ∓ and m s are the masses of the primary and the secondaryof inner and outer binaries. The values are shown in the scale-free unit with the gravitational constant, G = 1. m p m s a e θ φ ψ E P in 0.900 0.100 0.00100 0.900 1.50 0.00 0.00 3.14 1 . × − out 2.00 1.00 1.00 0.990 0.10 0.00 0.00 3.14 3.63 ( t ) H s d Figure 7.
The evolution of κ and H sd for the B-S system byusing the SD method. ( H , t ) ×10 ( H S D , t ) ×10 ( S D , t ) ×10 Figure 8.
The evolution of ε ( H ), ε ( H sd ) and ε (Γ sd ) for the B-Ssystem. In the ORG method, ε ( H ) and ε (Γ) are equivalent to ε ( H sd ) and ε (Γ sd ) with κ = 1 separately. To be convenient, weuse ε ( H sd ) and ε (Γ sd ) as the y -axis labels to represent both theSD and ORG cases (a similar style is used in other plots). Thedefinition of colors are the same as in Fig. 5. Notice in the plotof ε (Γ sd ), the scale of y -axis is two order of magnitude smallercompared to the upper plots. ( L x ) ×10 ( L y ) ×10 ( L z ) ×10 ( L ) ×10 Figure 9.
The cumulative error of three components of the totalangular momentum of the B-S system normalized to L (0) (initialvalue). The definition of colors are the same as in Fig. 5. and Γ sd (SD) within one integration step. Thus, ε (Γ sd ) rep-resents the numerical errors of the integrator.The angular momentum conservation is also checked inFig. 9. The normalized cumulative error of three compo-nents, (cid:101) ε ( L i ) = L i ( t ) − L i (0) | L (0) | ( i = x, y, z ) , (54)and the total one, (cid:101) ε ( L ) = (cid:113) [ L ( t ) − L (0)] | L (0) | , (55)are shown. Both the SD and the ORG methods provide asimilar level of relative numerical errors (10 − ) in the threecomponents of L and also in | L | . Besides, the error is inde-pendent of the variation of κ . This is consistent as we provedin Section 2.5.2. The mean motion resonance cannot be properly reproducedby the slow-down method. When there are multiple binaries,the orbital resonance between inner binaries can occur if c (cid:13)000
The cumulative error of three components of the totalangular momentum of the B-S system normalized to L (0) (initialvalue). The definition of colors are the same as in Fig. 5. and Γ sd (SD) within one integration step. Thus, ε (Γ sd ) rep-resents the numerical errors of the integrator.The angular momentum conservation is also checked inFig. 9. The normalized cumulative error of three compo-nents, (cid:101) ε ( L i ) = L i ( t ) − L i (0) | L (0) | ( i = x, y, z ) , (54)and the total one, (cid:101) ε ( L ) = (cid:113) [ L ( t ) − L (0)] | L (0) | , (55)are shown. Both the SD and the ORG methods provide asimilar level of relative numerical errors (10 − ) in the threecomponents of L and also in | L | . Besides, the error is inde-pendent of the variation of κ . This is consistent as we provedin Section 2.5.2. The mean motion resonance cannot be properly reproducedby the slow-down method. When there are multiple binaries,the orbital resonance between inner binaries can occur if c (cid:13)000 , 1–14 DAR integrator their periods are in resonance ratios. We investigate thiseffect by simulating a quadruple system (B-B). The initialcondition (Table 2) is generated by splitting the third bodyof the B-S system to a binary, which has the same periodas another. Here after the suffixes “in1” and “in2” indicatethe two inner binaries. t KL of the two binaries are about 45and 88. We evolve the system to t = 160 with the same ds used in the integration of the B-S system. The ORG methodtakes about 8 . × steps while the SD method takes about9 . × steps (9 times faster).The orbital evolution is shown in Fig. 10. t KL , in2 ofthe simulated result is a much larger than the predictionof Eq. 53. The behaviour of the first inner binary is sim-ilar to the case of the B-S system. Compared to the caseof B-S system, the cumulative divergence of a for the SDand ORG methods becomes obvious. After about two t KL ,the relative difference is the order of 10 − for a in1 and a in2 .However, this error can be neglected since the Kozai-Lidoveffect dominates the secular evolution. Thus, the SD methodcan still provide a reasonable result.On the other hand, since a of two inner binaries donot evolve significantly, it is easy to observe the differencecaused by the orbital resonance between the SD and theORG method. Although the SD method cannot keep theoriginal period ratio of the two binaries, the mean motionresonance still exist because κ i of the two binaries are relatedby Eq. 22. The upper panel of Fig. 11 shows the variationof κ i , where κ in2 is twice of κ in1 . Therefore, the period ratiochanges from 1 : 1 to 1 : 2. Since the orbits are slowed down,the resonance becomes stronger, which results in the largerdrift in the SD case. ε (Γ sd ) keeps the order of 10 − . Both SD and ORG meth-ods show jumps when e in1 passes the maximum value due tothe large numerical error at the pericenter of the first innerbinary. Such error can be reduced with a smaller ∆ s (notshown here). The behaviour of | L | is independent of e in1 andits relative error has an order of 10 − . In star clusters, hyperbolic encounters between stars andbinaries 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 aperturber is far away, the binary has a weak perturbation sothat κ is large. Once the encounter starts, κ drops fast. Thus,it is necessary to limit the change rate of κ , for example, byusing the timescale criterion (Eq. 28).We test a hyperbolic encounter between a massive bi-nary and a low-mass binary (HB-B). The initial conditionis shown in Table 3. The mass ratio between the two innerbinaries is 100. The hyperbolic (outer) orbit has an initialseparation of 2 .
86 and a pericenter distance of 0 . a in1 ). Initially E out = − .
00 so that the time toreach the pericenter is about 1 .
71. Based on Eq. 22, κ in1 hasthe maximum value of about 3 × and can reaches theminimum limit of 1 .
0. Thus, κ in1 varies about four order ofmagnitude in one encounter. By using the timescale crite-rion with c = 0 . κ in1 is limited to about 10 ,which is one magnitude smaller.Due to the high mass ratio, the second inner binaryfeel a strong perturbation from the first inner binary. The final fate of the second binary depends on its orbital phaseduring the encounter. In the ORG case, three values of E in2 are used to test the effect of orbital phases. The SD modeluse the third value of E in2 listed in Table 3.The result of orbital revolution is shown in Fig. 12. Af-ter the encounter, a jump appears in all orbital parameters.The initial E in2 significantly influences the final orbit of thesecondary inner binary (a large divergence appears). How-ever, the SD and ORG method (SD-3 and ORG-3) can pro-vide a similar result when the initial E in2 are same (3 . κ in1 (for the SD-3model) calculated by the perturbation criterion (Eq. 22) andby both the perturbation and the timescale criterion. Thetimescale criterion ensures that κ in1 changes more smoothlyduring the encounter. ε ( H sd ) has a larger error of (10 − )compared to ε ( H ) at the beginning because κ in1 is large.It would be worse if the timescale criterion is not applied.Although ε (Γ sd ) is not small, | L | is still well conserved (er-ror has an order of 10 − ), which again suggests that theconservation of L is independent of κ .With ∆ s ≈ . × − , the ORG method requiresabout 2 . × steps to reach t = 4 while the SD methodonly needs 3 . × steps. Thus, the SD method provides a60 times faster performance. Notice here we only evolve thesystem in a short time interval around the time of encounter.In a star cluster, most of the time the binaries are weaklyperturbed, thus averagely κ i (cid:29) . Therefore, the totalintegration steps of binaries are significantly reduced dur-ing the long-term evolution. This example suggests that inthe simulation of star clusters with many binaries, the SDmethod can dramatically improves the performance whilethe statistical result of encounters can be well reproduced. In this work, we mathematically and numerically describethe slow-down time-transformed explicit symplectic (SDAR)method. It combines the benefit of the symplectic integra-tor which conserves the Hamiltonian and angular momen-tum and the high efficiency of the slow-down method tohandle the long-term evolution of hierarchical systems andclose encounters. An implementation of the method writtenin the c++ language, sdar , is publicly available. The codemodularized for easily plugging in other N -body codes.The Hamiltonian and the corresponding equation ofmotion are discussed in detail. Although the physical energy( H ) is not conserved, the method can provide a correct sec-ular evolution. We mathematically prove that L is alwaysconserved with the slow-down method. We also discussedhow to measure the numerical error, ε (Γ sd ), in the absenceof energy conservation.We show that in the LogH method, for a Kepler orbitthe integration step d s has the geometric meaning of eccen-tric anomaly(Eq. 47). Using this feature, we can determinethe number of steps per orbit by setting ds using Eq. 47.This is very useful for controlling the integration error. c (cid:13) , 1–14 Long Wang et al.
Table 2.
The initial condition of the hierarchical quadruple (B-B) system. m p m s a e θ φ ψ E P in1 0.900 0.100 0.00100 0.900 1.50 0.00 0.00 3.14 1 . × − in2 1.800 0.200 0.00126 0.900 0.10 0.00 0.00 3.14 1 . × − out 2.00 1.00 1.00 0.990 0.100 0.00 0.00 3.14 3.63 a ×10 in1
10 ×10 +1.261×10 in2 out e Figure 10.
The evolution of orbital parameters of the two inner and outer binaries for the B-B system. The plotting style is similar toFig. 5. ( t ) ×10 in1in21.51.00.50.0 ( S D , t ) ×10 SDORG0 25 50 75 100 125 150Time024 ( L ) ×10 Figure 11.
The evolution of κ , ε (Γ sd ) and (cid:101) ε ( L ) for the B-Bsystem. The upper panel shows κ i of the two inner binaries. Themiddle and lower panels compared ε (Γ sd ) and (cid:101) ε ( L ) for the SDand ORG cases. On the other hand, when the LogH method is imple-mented in a hybrid N -body integrator, it is necessary toensure that the integration can stop at a given time, in or-der to calculate the interactions between the members offew-body systems and the global system. However, becausetime is an integrated value, before the integration finishes,the next time is unknown. A few iteration of integration isneeded to converge the next time to a certain value with agiven limit of error. Such iteration can be expensive if thetime synchronization is frequent. Thus, in order to avoidtoo many synchronization requests, it is necessary to de-sign a good criterion for determination of the subsystemsthat need to be handled by the SDAR method. Especially,the number of iteration steps should be much less than thenumber of integration steps.By using the ar code, we show three numerical ex-amples (B-S, B-B and HB-B) to indicate that the SDARmethod can well reproduce the Kozai-Lidov oscillation andcan give an acceptable result of a hyperbolic encounter be-tween two binaries. When the averaged value of κ is large,the slow-down method can provide significant performanceimprovement. Especially, for binaries with weak perturba-tions, the method can provide a few order of magnitude lessintegration steps. Thus, by combining this algorithm in anhybrid integrator for simulating a large N -body system, itis expected that a large fraction of binaries and hierarchi-cal systems can be efficiently handled. Such a code will beavailable soon in our following-up work. c (cid:13)000
The evolution of κ , ε (Γ sd ) and (cid:101) ε ( L ) for the B-Bsystem. The upper panel shows κ i of the two inner binaries. Themiddle and lower panels compared ε (Γ sd ) and (cid:101) ε ( L ) for the SDand ORG cases. On the other hand, when the LogH method is imple-mented in a hybrid N -body integrator, it is necessary toensure that the integration can stop at a given time, in or-der to calculate the interactions between the members offew-body systems and the global system. However, becausetime is an integrated value, before the integration finishes,the next time is unknown. A few iteration of integration isneeded to converge the next time to a certain value with agiven limit of error. Such iteration can be expensive if thetime synchronization is frequent. Thus, in order to avoidtoo many synchronization requests, it is necessary to de-sign a good criterion for determination of the subsystemsthat need to be handled by the SDAR method. Especially,the number of iteration steps should be much less than thenumber of integration steps.By using the ar code, we show three numerical ex-amples (B-S, B-B and HB-B) to indicate that the SDARmethod can well reproduce the Kozai-Lidov oscillation andcan give an acceptable result of a hyperbolic encounter be-tween two binaries. When the averaged value of κ is large,the slow-down method can provide significant performanceimprovement. Especially, for binaries with weak perturba-tions, the method can provide a few order of magnitude lessintegration steps. Thus, by combining this algorithm in anhybrid integrator for simulating a large N -body system, itis expected that a large fraction of binaries and hierarchi-cal systems can be efficiently handled. Such a code will beavailable soon in our following-up work. c (cid:13)000 , 1–14 DAR integrator 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 toreach the pericenter. The inner binary orbits are the same as in the B-B system except in the ORG models, E in2 uses three values fortesting the impact of initial orbital phases. The SD model (SD-3) chooses the third value of E in2 . m p m s a e θ φ ψ E P in1 0.900 0.100 0.00100 0.900 1.50 0.00 0.00 3.14 1 . × − in2 0.00900 0.00100 0.00200 0.900 1.50 0.00 0.00 3.00, 3.50, 3.14 5 . × − out 1.00 0.0100 -1.00 1.02500 0.100 0.00 0.00 -2.00 1.71 a ×10 +10 in1 ORG-1 ORG-2 ORG-3 SD-305 ×10 in2 out e ×10 +9×10
05 ×10 +3.1415 2.50.02.5 50 ×10
10 ×10
05 ×10
02 ×10 +1.5 0 1 2 3 4Time0.51.0 0 1 2 3 4Time1.0001.005 ×10 Figure 12.
The evolution of orbital parameters of the two inner binaries and the outer encounter for the HB-B system. The solid anddashed lines represent the ORG and SD methods separately. Colors indicate different E in2 in the ORG cases. The SD-3 and ORG-3 havethe same initial E in2 . The plotting style is similar to Fig. 5. ACKNOWLEDGMENTS
L.W. thanks the financial support from JSPS Interna-tional Research Fellow (School of Science, The universityof Tokyo).
REFERENCES
Aarseth S. J., 2003, Gravitational N-Body Simulations,Cambridge University PressAntognini J. M. O., 2015, MNRAS, 452, 3610Antonini F., Chatterjee S., Rodriguez C. L., Morscher M.,Pattabiraman B., Kalogera V., Rasio F. A., 2016, ApJ,816, 65Askar A., Szkudlarek M., Gondek-Rosi´nska D., Giersz M.,Bulik T., 2017, MNRAS, 464, L36Bailyn C. D., 1995, ARA&A, 33, 133Banerjee S., Baumgardt H., Kroupa P., 2010, MNRAS, 402,371Breen P. G., Heggie D. C., 2013, MNRAS, 432, 2779Binney J., Tremaine S., 1987, Princeton, NJ, PrincetonUniversity PressBurdet C. A., 1967, ZaMP, 18, 434Burdet C. A., 1968, ZaMP, 19, 345Davies M. B., Piotto G., de Angeli F., 2004, MNRAS, 349,129 Di Carlo U. N., Giacobbo N., Mapelli M., Pasquato M.,Spera M., Wang L., Haardt F., 2019, MNRAS, 487, 2947Fujii M. S., Portegies Zwart S., 2011, Sci, 334, 1380Gualandris A., Portegies Zwart S., Sipior M. S., 2005, MN-RAS, 363, 223Gvaramadze V. V., Gualandris A., Portegies Zwart S.,2009, MNRAS, 396, 570Hairer, E., 1997, Applied Numerical Mathematics, 25, 219-227Harrington R. S., 1968, AJ, 73, 190Heggie D. C., 1973, ASSL, 34, ASSL...39Heggie D. C., 1975, MNRAS, 173, 729Heggie D. C., Giersz M., 2008, MNRAS, 389, 1858H´enon M., 1961, AnAp, 24, 369Hida Y., Li X. S., Bailey D. H., 2001, ”Algorithms for quad-double precision floating point arithmetic,” Proceedings15th IEEE Symposium on Computer Arithmetic. ARITH-15 2001, Vail, CO, USA, 2001, pp. 155-162.Hills J. G., 1975, AJ, 80, 809Hurley J. R., Pols O. R., Aarseth S. J., Tout C. A., 2005,MNRAS, 363, 293Hypki A., Giersz M., 2013, MNRAS, 429, 1221Kozai Y., 1962, AJ, 67, 591Kustaanheimo P., Stiefel E., 1965, J. Reine Angew. Math.,218, 204Leigh N., Sills A., Knigge C., 2011, MNRAS, 416, 1410 c (cid:13) , 1–14 Long Wang et al. ( t ) ×10 in1 (used) in1 (pert)1.00.50.00.51.0 ( S D , t ) ×10 SD-3ORG-30 1 2 3 4Time024 ( L ) ×10 Figure 13.
The evolution of κ in1 (SD-3), ε (Γ sd ) and (cid:101) ε ( L ) forthe HB-B system. In the upper panel, the dashed curve (used) iscalculated by both the perturbation and the timescale criterionand the solid curve (pert) only applies the perturbation criterion.The SD-3 and ORG-3 models are compared for ε (Γ sd ) and (cid:101) ε ( L )with a similar plotting style to that of Fig. 11. Leonard P. J. T., Duncan M. J., 1988, AJ, 96, 222Leonard P. J. T., Duncan M. J., 1990, AJ, 99, 608Lidov M. L., 1962, P&SS, 9, 719Mikkola S., Aarseth S. J., 1996, CeMDA, 64, 197Mikkola S., Tanikawa K., 1999, MNRAS, 310, 745Mikkola S., Aarseth S., 2002, CeMDA, 84, 343Naoz S., Farr W. M., Lithwick Y., Rasio F. A., TeyssandierJ., 2013, MNRAS, 431, 2155Naoz S., 2016, ARA&A, 54, 441Spitzer L., 1987, Dynamical evolution of globular clusters,Princeton University PressO’Leary R. M., Rasio F. A., Fregeau J. M., Ivanova N.,O’Shaughnessy R., 2006, ApJ, 637, 937Portegies Zwart S. F., McMillan S. L. W., 2000, ApJL, 528,L17Poveda A., Ruiz J., Allen C., 1967, BOTT, 4, 86Preto M., Tremaine S., 1999, AJ, 118, 2532Wang L., 2020, MNRAS, 491, 2413Yoshida H., 1990, PhLA, 150, 262Yu Q., Tremaine S., 2003, ApJ, 599, 1129
APPENDIX A: PROOF FOR THECONSERVATION OF ANGULAR MOMENTUM
Here we provide the proof for the conservation of angularmomentum in the slow-down method. First, the conserva-tion of angular momentum for an isolated binary can be described as { L b , H b } = L b = m r × v + m r × v H b = 12 m v + 12 m v − Gm m | r − r | (A1)where L b , H b are the angular momentum and Hamiltonianof a binary system in the center-of-the-mass frame.After a slow-down factor κ ( t ) is included, { L b , H sd } = { L b , κ H b } = 1 κ { L b , H b } = . (A2) L b is still conserved. When the center-of-the-mass term isincluded, H sd = 1 κ (cid:20) m ( v − v cm ) + 12 m ( v − v cm ) − Gm m | r − r | (cid:21) + 12 ( m + m ) v , (A3)it can be shown that the conservation (Eq. A1) still exists.If a new body is added to form a triple, the additionalterm appears in the angular momentum: L = L b + L . (A4)Then { L , H b } = { L b , H b } + { L , H b } = { L , H b } . (A5)Since L is independent of binary position and velocity, { L , H b } = { L , H b } = . (A6)The additional term to H sd from the third body is H = 12 m v − (cid:88) i Gm i m | r i − r | . (A7)It can be proved that { L , H } = . (A8)Thus, { L , H sd } = { L , H b } + { L , H } = . (A9)The angular momentum is conserved for the triple.For systems with one binary and many singles, 4 th , th ... particles can be added one by one, and for the k th particle,the additional Hamiltonian term is H i = 12 m i v i − i − (cid:88) j Gm i m j | r i − r j | . (A10)Since the second term in H i is a linear combination of k − { L , H i } = . (A11)Eventually, { L , H sd } = (A12)Thus, L is conserved in such system.On the other hand, here we add many singles to a slow-down binary and indicate that L is conserved. We can con-sider this process in an opposite way: add one slow-downbinary to a systems of many singles, the new L is still con-served. Thus, it can be proved that adding arbitrary slow-down binaries to the systems, L is always conserved. c (cid:13)000