IMPLEMENTING FEW-BODY ALGORITHMIC REGULARIZATION WITH POST-NEWTONIAN TERMS

After the introduction of electronic computers, those who carried out simulations of the gravitational N-body problem soon realized that the classical methods of numerical integration were often not satisfactorily accurate due simply to the strong 1/r² character of the gravitational force.

The situation changed when Kustaanheimo & Stiefel (1965) published their Kustaanheimo–Stiefel (KS) transformation from four-dimensional to three-dimensional space. The special case of a planar system had been known for a long time (Levi-Civita 1920) but it had turned out that a similar transformation in the three-dimensional space was not possible. The situation for the general N-body problem improved only after the KS transformation became well known, due in large part to publication by Stiefel & Scheifele (1971) of their text, which comprehensively discussed the application of the KS transformation to the perturbed two-body problem. Later, Aarseth & Zare (1974), Heggie (1974), Zare (1974), and Mikkola & Aarseth (1993) applied the KS transformation to the general three-body and N-body problems.

An entirely new way of regularizing close encounters was invented simultaneously by Mikkola & Tanikawa (1999a, 1999b) and by Preto & Tremaine (1999). This new method introduced the so-called logarithmic Hamiltonian (LogH). Together with the simple leapfrog algorithm, this new method gives regular results for close encounters; in fact a correct trajectory is obtained for the two-body problem. Other similar methods are described by Huang & Leimkuhler (1997).

The remaining problem was that none of these regularization methods could be easily applied to systems with extremely large mass ratios. In an attempt to solve this problem, Mikkola & Aarseth (2002) introduced the time-transformed leapfrog (TTL). This method is in some cases mathematically equivalent to the LogH method, but it is more general, and arbitrary mass ratios are allowed. The drawback of this method, however, is that in some cases the roundoff error affects the results considerably. In these new methods (LogH, TTL) the regularization is achieved by using the leapfrog, hence the name "algorithmic regularization." More details about the thus-far mentioned methods can be found in the book by Aarseth (2003).

Since the leapfrog alone is rarely accurate enough, one must supplement the method with the extrapolation algorithm (Gragg 1964, 1965; Bulirsch & Stoer 1966; Press et al. 1986), or with higher-order leapfrogs (Yoshida 1990), in order to get highly accurate results. The efficiency of the extrapolation procedure requires that the basic algorithm (leapfrog in algorithmic regularization, modified midpoint method in the KS-regularized codes) have a certain symmetry. In the case of the leapfrog this means time reversibility. If the system has velocity-dependent forces, such as relativistic post-Newtonian (PN) terms (Soffel 1989), then the required symmetry is more difficult to obtain. One way to solve this problem was recently obtained by Mikkola & Merritt (2006), who formally doubled the dimensionality of the parameter space and constructed a generalized midpoint method (GAR). This algorithm allows the use of algorithmic regularization and velocity-dependent perturbations. However, these authors discussed the new algorithm essentially only for the case of the perturbed two-body problem.

In this paper, we discuss the details of our most recent implementation of the algorithmic regularization method. This uses the chain structure, the same structure as in the KS-CHAIN algorithm of Mikkola & Aarseth (1993). That device, together with a new time-transformation function, significantly reduces the roundoff problems and makes the new code a good alternative for simulations of strongly-interacting few-body systems. In the final section, we present some applications of the chain routine to the problem of relativistic orbits around the supermassive black hole (SMBH) at the Galactic center.

In this paper, we use the following basic symbols:

$$\begin{equation*} \begin{array}{@{}l@{\qquad}l@{}} G=1 & {\rm gravitational\ constant}\\ \ms t & {\rm time}\\ \ms m_k & {\rm mass\ of\ body}\ k \\ \ms {\fat p}_k & {\rm momentum} \\ \ms \fat r_k & {\rm position} \\ \ms \fat v_k=\dot{\fat r_k} & {\rm velocity}\\ \ms T=\sum _{k=1}^N{\fat p_k^2}\big/{2m_k} & {\rm kinetic\ energy}\\ \ms U=\sum _{0<i<j\le N}{m_im_j}\!/{|\fat r_i-\fat r_j|} & {\rm potential\ energy}\\ \ms E=T-U & {\rm total\ energy}\\ \ms B=-E=U-T & {\rm binding\ energy} \end{array} \end{equation*}$$

In the few-body problem, close approaches of two bodies are common, and these events require highly accurate computations since the energies involved are large. Thus, a method that can accurately advance the motions of few-body systems must be accurate and efficient for the two-body problem. In addition to the classical KS transformation there are some new algorithms that satisfy this requirement.

An additional, and most important, problem is the roundoff error. This becomes a major problem if, e.g., center-of-mass coordinates are used and there are close binaries and/or close encounters in the simulated system. The solution to this problem is the use of the chain structure, originally applied with the KS transformation (Mikkola & Aarseth 1993) in order to KS-regularize all the short(est) distances. Later it became clear that the chain structure was also beneficial in significantly reducing the roundoff error.

In this section, we first review the ingredients we need in the construction of the final N-body algorithm. These basic algorithms produce regular results in close approaches and their results can be improved to high precision using, e.g., the Bulirsch & Stoer (1966) extrapolation method.

3.1. LogH, TTL, and GAR

3.1.1. LogH

Recently, Mikkola & Tanikawa (1999a) and Preto & Tremaine (1999) pointed out that the LogH in extended phase space,

$$\begin{equation} {\Lambda }=\ln (T+B)-\ln (U), \end{equation} \tag{ 1 }$$

where B (the binding energy) is the momentum of time, gives the equations of motion in the form

$$\begin{equation} t^{\prime }={\frac{\partial {\Lambda }}{\partial B}}=1/(T+B) \end{equation} \tag{ 2 }$$

$$\begin{equation} \fat r_k^{\prime }={\frac{\partial T}{\partial \fat p_k}}\bigg/(T+B) \end{equation} \tag{ 3 }$$

$$\begin{equation} B^{\prime }={\frac{\partial U}{\partial t}}\bigg/U \end{equation} \tag{ 4 }$$

$$\begin{equation} \fat p_k^{\prime}=\doo{U}{\fat r_k}\bigg/U. \end{equation} \tag{ 5 }$$

Here we include Equation (4) although this is needed only if there is a time-dependent potential to be added to the N-body potential U. Since the right-hand sides of these equations do not depend on the left-hand-side variables, a leapfrog algorithm is possible. This may be symbolized as

$$\begin{equation} {\fat X}(h/2){\fat V}(h){\fat X}(h)\cdots {\fat V}(h){\fat X}(h/2), \end{equation} \tag{ 6 }$$

where ${\fat X}(s)$ means solution for the coordinate Equations (2), (3) with constant B and $\fat p_k$ over an integration step of length = s:

$$\begin{equation} \delta t=s/(T+B);\quad t\rightarrow t+\delta t;\quad \fat r_k\rightarrow \fat r_k+\delta t \: {\frac{\partial T}{\partial \fat p_k}}. \end{equation} \tag{ 7 }$$

Correspondingly ${\fat V}(s)$ signifies the operation

$$\begin{equation} \widetilde{\delta t}=s/U;\quad B\rightarrow B+\widetilde{\delta t} \: {\frac{\partial U}{\partial t}}; \quad \fat p_k\rightarrow \fat p_k +\widetilde{\delta t} {\frac{\partial U}{\partial \fat r_k}}, \end{equation} \tag{ 8 }$$

which solves Equations (4) and (5) for constant t and $\fat r_k$.

For the case of two bodies only, this algorithm produces correct trajectories, with only an O(h³) phase error. This is true even for the collision orbit and the energy conservation is, in this case, typically at the machine precision level.

3.1.2. TTL

The second important ingredient in our algorithm is the TTL (Mikkola & Aarseth 2002). The basic idea is to introduce a coordinate-dependent time-transformation function $\Omega (\ldots,\fat r_k,\ldots)$ and a new variable ω which is supposed to have the same numerical value as Ω, but that value is obtained from the differential equation

$$\begin{equation} \dot{\omega }=\dot{\Omega }=\sum _k{\frac{\partial \Omega }{\partial \fat r_k}}\cdot { \fat v_k}. \end{equation} \tag{ 9 }$$

The equations of motion can be written as

$$\begin{equation} t^{\prime }=1/\omega \end{equation} \tag{ 10 }$$

$$\begin{equation} \fat r_k^{\prime }=\fat v_k/\omega \end{equation} \tag{ 11 }$$

$$\begin{equation} \fat v_k^{\prime }=\fat A_k/\Omega \end{equation} \tag{ 12 }$$

$$\begin{equation} \omega ^{\prime }=\sum _k{\frac{\partial \Omega }{\partial \fat r_k}}\cdot { \fat v_k}/\Omega, \end{equation} \tag{ 13 }$$

where $\fat A_k$ is the acceleration of the kth particle. The structure of these equations allows the construction of a leapfrog algorithm:

$$\begin{eqnarray} &&\hskip-15pc{{\fat X(s)}:}\ \ \ \ \ \nonumber\\ &&\delta t=s/\omega; \quad t\rightarrow t+\delta t;\quad \fat r_k\rightarrow \fat r_k+\delta t \: \fat v_k \end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray} &&\hskip-11pt{{ \fat V(s)}:} \nonumber\\ &&\quad\widetilde{\delta t}=s/\Omega; \quad \fat v_k\rightarrow \fat v_k +\widetilde{\delta t} \fat A_k; \quad \omega \rightarrow \omega +\widetilde{\delta t} \: \langle \dot{\Omega }\rangle,\nonumber\\ \end{eqnarray} \tag{ 15 }$$

where $\langle \dot{\Omega }\rangle$ is the average of this quantity over the step, i.e.,

$$\begin{equation} \langle \dot{\Omega }\rangle =\sum _k{\frac{\partial \Omega }{\partial \fat r_k}}\cdot \big(\fat v_k^{\rm old}+\fat v_k^{\rm new}\big)\big/2. \end{equation} \tag{ 16 }$$

Here, the superscripts "old" and "new" refer to $\fat v_k$ values before and after the velocity advancement in the operation (15).

3.1.3. GAR

Here we concisely review the GAR method. Following Mikkola & Merritt (2006), we consider a differential equation

$$\begin{equation} \dot{\fat z}=\fat f(\fat z), \end{equation} \tag{ 17 }$$

and an approximation for its solution (over a short step = h) written as

$$\begin{equation} \fat z(h)\approx \fat z(0)+\fat d(\fat z(0),\,h). \end{equation} \tag{ 18 }$$

Here the increment $\fat d(\fat z_0,h)$ can be any approximation suitable for the particular equation in question. If one writes, instead of Equation (17), the two equations

$$\begin{equation} \dot{\fat x}=\fat f(\fat y);\qquad \dot{\fat y}=\fat f(\fat x), \end{equation} \tag{ 19 }$$

and solves this pair with the initial values $\fat x(0)=\fat y(0)=\fat z(0)$, the solution obviously is $\fat x(t)=\fat y(t)=\fat z(t)$. Using the pair of equations one can write a leapfrog as

$$\begin{equation*} \nonumber \fat x_{\scriptstyle {\frac{1}{2}}}=\fat x_0+\frac{h}{2}\fat f(\fat y_0); \quad \fat y_1=\fat y_{\scriptstyle {\frac{1}{2}}}+h\fat f\big(\fat x_{\scriptstyle {\frac{1}{2}}}\big); \quad \fat x_1=\fat x_{\scriptstyle {\frac{1}{2}}}+\frac{h}{2}\fat f(y_1), \end{equation*}$$

which is actually nothing but the well-known modified midpoint method. In the above, one can split the advancement of y into two operations to get

$$\begin{equation} \fat x_{\scriptstyle {\frac{1}{2}}}=\fat x_0+\frac{h}{2}\fat f(\fat y_0); \qquad \fat y_{\scriptstyle {\frac{1}{2}}}=\fat y_0+\frac{h}{2} \fat f\big(\fat x_{\scriptstyle {\frac{1}{2}}}\big) \end{equation} \tag{ 20 }$$

$$\begin{equation} \fat y_1=\fat y_{\scriptstyle {\frac{1}{2}}}+\frac{h}{2}\fat f\big(\fat x_{\scriptstyle {\frac{1}{2}}}\big); \qquad \fat x_1=\fat x_{\scriptstyle {\frac{1}{2}}}+\frac{h}{2}f(y_1), \end{equation} \tag{ 21 }$$

and this can be readily generalized using the more general increment $\fat d(\fat z,\,h)$. This results in the generalized midpoint method

$$\begin{equation} \fat x_{\scriptstyle {\frac{1}{2}}}=\fat x_0+\fat d\left(\fat y_0,\,\frac{h}{2}\right); \qquad \fat y_{\scriptstyle {\frac{1}{2}}}=\fat y_0-\fat d\left(\fat x_{\scriptstyle {\frac{1}{2}}},\,{-}\frac{h}{2}\right), \end{equation} \tag{ 22 }$$

$$\begin{equation} \fat y_1=\fat y_{\scriptstyle {\frac{1}{2}}}+\fat d\left(\fat x_{\scriptstyle {\frac{1}{2}}},\,\frac{h}{2}\right); \qquad \fat x_1=\fat x_{\scriptstyle {\frac{1}{2}}}-\fat d\left(\fat y_1,\,{-}\frac{h}{2}\right). \end{equation} \tag{ 23 }$$

This method has the great advantage that one can use any special approximation $\fat d$ and the algorithm is time reversible (albeit only in the extended $\fat x\fat y$-space) and thus suitable to be used as the basic integrator in an extrapolation method. Specifically, one may stress that the increment $\fat d$ may be computed using LogH or TTL which produce good approximations in the case of close encounters. In addition, the appearance of velocity-dependent forces is not problematic here, as shown by Mikkola & Merritt (2006).

While in the algorithms discussed above one can use any coordinate system, the roundoff error can be a major problem in the case of close encounters if the coordinates of the approaching bodies are measured from a distant origin. This problem is significantly reduced by utilizing the chain structure. This was originally used in a KS-regularization algorithm (Mikkola & Aarseth 1993) in order to get all the short distances regularized by the KS transformation, but here the sole purpose is the roundoff reduction, for which the device has proved itself.

4.1. Basic Formulation

Let $T={1\over 2}\!\sum _{k=1}^N m_k\fat v_k^2$ be the kinetic energy and U = ∑_i<j⩽Nm_im_j/r_ij the potential such that the total energy is E = T − U; the binding energy is B = U − T.

One forms a chain of particles such that the shortest relative vectors are in the chain (Mikkola & Aarseth 1993). We stress again that the main purpose of using the chain structure in this method is to reduce the (often significant) effects of the roundoff error.

Let us collect the chain coordinates $\fat X_k=\fat r_{i_k} -\fat r_{j_k}$ in the vector $\fat X=(\fat X_1,\,\fat X_2,\ldots,\fat X_{N-1})$ and let the corresponding velocities $\fat V_k=\fat v_{i_k}-\fat v_{j_k}$ be in the vector $\fat V=(\fat V_1,\fat V_2,\ldots,\fat V_{N-1})$. Then, the Newtonian equations of motion may be formally written as

$$\begin{equation} \dot{\fat X}=\fat V \end{equation} \tag{ 24 }$$

$$\begin{equation} \dot{\fat V}=\fat A(\fat X)+\fat f, \end{equation} \tag{ 25 }$$

where $\fat A$ is the N-body acceleration and $\fat f$ is some external acceleration (e.g. due to other bodies).

One may use the two equivalent time transformations (Mikkola & Merritt 2006)

$$\begin{equation} ds=[\alpha (T+B)+\beta \omega +\gamma ]\,dt = [\alpha U+\beta \Omega +\gamma ]\,dt, \end{equation} \tag{ 26 }$$

where s is the new independent variable, α, β, and γ are adjustable constants, Ω is an optional function of the coordinates $\Omega =\Omega (\fat X)$, while the initial value ω(0) = Ω(0) and the differential equation

$$\begin{equation} \dot{\omega }={\frac{\partial \Omega }{\partial \fat X}}\cdot \fat V \end{equation} \tag{ 27 }$$

determine the values of ω (in fact, ω(t) = Ω(t) along the exact solution).

The time transformation thus introduced regularizes the two-body collisions if one uses the simple leapfrog algorithm as a basic integrator, results from which can, and must, be improved using an extrapolation method (e.g. Bulirsch & Stoer 1966; Press et al. 1986).

It is possible to divide the equations of motion into two categories (where derivatives with respect to the new independent variable s are denoted by a prime):

Coordinate equations:

$$\begin{equation} t^{\prime }=1/(\alpha (T+B)+\beta \omega +\gamma) \end{equation} \tag{ 28 }$$

$$\begin{equation} \fat X^{\prime }=t^{\prime }\fat V. \end{equation} \tag{ 29 }$$

Velocity equations:

$$\begin{equation} \tilde{t}^{\prime }=1/(\alpha U+\beta \Omega +\gamma) \end{equation} \tag{ 30 }$$

$$\begin{equation} \fat V^{\prime }=\tilde{t}^{\prime }(\fat A+\fat f) \end{equation} \tag{ 31 }$$

$$\begin{equation} \omega ^{\prime }=\tilde{t}^{\prime }{\frac{\partial \Omega }{\partial \fat X}}\cdot \fat V \end{equation} \tag{ 32 }$$

$$\begin{equation} B^{\prime }=-\tilde{t}^{\prime }{\frac{\partial T}{\partial \fat V}}\cdot \fat f. \end{equation} \tag{ 33 }$$

In these equations the right-hand sides do not depend on the variables on the left. Consequently it is possible to construct a regular leapfrog algorithm for obtaining the solutions. The leapfrog results then can easily be improved with the extrapolation method.

4.2. Details

4.2.1. Finding and Updating the Chain

First we find the shortest interparticle vector which is adopted as the first part of the chain. The chain is then augmented by adding the relative vector to the particle nearest to one or the other end of the existing chain. When all particles are included, they are re-numbered along the chain as 1, 2, ..., N for ease of programming.

To reduce roundoff problems, the transformation from the old chain vectors ${\fat X}_k$ to the new ones is done directly by expressing the new chain vectors as sums of the old ones as in Mikkola & Aarseth (1993).

4.2.2. Transformations

When the particles are renamed along the chain as 1, 2, ..., N one can evaluate

$$\begin{equation} \fat X_k=\fat r_{k+1}-\fat r_k \end{equation} \tag{ 34 }$$

$$\begin{equation} \fat V_k=\fat v_{k+1}-\fat v_k. \end{equation} \tag{ 35 }$$

The center-of-mass quantities are

$$\begin{equation} M=\sum _k m_k \end{equation} \tag{ 36 }$$

$$\begin{equation} \fat r_{\rm cm}=\sum _k m_k \fat r_k/M \end{equation} \tag{ 37 }$$

$$\begin{equation} \fat v_{\rm cm}=\sum _k m_k \fat v_k/M. \end{equation} \tag{ 38 }$$

The inverse transformation is done by simple summation

$$\begin{equation} \widetilde{\fat r}_1 ={\bf 0} \end{equation} \tag{ 39 }$$

$$\begin{equation} \widetilde{\fat v}_1 ={\bf 0} \end{equation} \tag{ 40 }$$

$$\begin{equation} \widetilde{\fat r}_{k+1} =\widetilde{\fat r}_k + \fat X_k \end{equation} \tag{ 41 }$$

$$\begin{equation} \widetilde{\fat v}_{k+1} =\widetilde{\fat v}_k + \fat V_k, \end{equation} \tag{ 42 }$$

followed by reduction to the center of mass

$$\begin{equation} \widetilde{\fat r}_{\rm cm}=\sum _k m_k \widetilde{\fat r}_k/M \end{equation} \tag{ 43 }$$

$$\begin{equation} \widetilde{\fat v}_{\rm cm}=\sum _k m_k \widetilde{\fat v}_k/M \end{equation} \tag{ 44 }$$

$$\begin{equation} \fat r_{k}=\widetilde{\fat r}_k-\widetilde{\fat r}_{\rm cm} \end{equation} \tag{ 45 }$$

$$\begin{equation} \fat v_{k}=\widetilde{\fat v}_k-\widetilde{\fat v}_{\rm cm}. \end{equation} \tag{ 46 }$$

Note that it is not always necessary to reduce the coordinates to the center-of-mass system since accelerations only depend on the differences.

4.2.3. Equations of Motion and the Leapfrog

One writes the equations of motion as

$$\begin{equation} \dot{\fat X}_k=\fat V_k \end{equation} \tag{ 47 }$$

$$\begin{equation} \dot{\fat V}_k=\fat F_{k+1}-\fat F_k+\fat f_{k+1}-\fat f_k, \end{equation} \tag{ 48 }$$

where $\fat f_k$ are the individual external accelerations and the N-body accelerations $\fat F_k$ are

$$\begin{equation} \fat F_k=-\sum _{j\ne k}m_j\frac{\fat r_{jk}}{|\fat r_{jk}|^3}, \end{equation} \tag{ 49 }$$

where, for j < k,

$$\begin{equation} \fat r_{jk}=\left\{ \begin{array}{@{}l@{\qquad}l@{}} \fat r_k-\fat r_j;& {\rm if}\quad k>j+2 \\ \fat X_j;& {\rm if}\quad k=j+1 \\ \fat X_j+\fat X_{j+1};& {\rm if}\quad k=j+2.\end{array}\right. \end{equation} \tag{ 50 }$$

For k > j one uses $\fat r_{jk}=-\fat r_{kj}$. In the acceleration computation the use of $\fat X_j$ and $\fat X_j+\fat X_{j+1}$ reduces the roundoff effect significantly. This is one of the most important features of the algorithm.

The kinetic energy is evaluated as usual

$$\begin{equation} T=\frac{1}{2} \sum _k m_k \fat v_k^2, \end{equation} \tag{ 51 }$$

while the potential

$$\begin{equation} U=\sum _{i<j}\frac{m_im_j}{|\fat r_{ij}|} \end{equation} \tag{ 52 }$$

is obtained along with the accelerations according to Equation (50). For the time transformation function Ω it seems advantageous to use

$$\begin{equation} \Omega =\sum _{i<j}\frac{\Omega _{ij}}{ |\fat r_{ij}|}. \end{equation} \tag{ 53 }$$

Here Ω_ij are adjustable constants (see below).

Now one is able to evaluate the two time transformation functions

$$\begin{equation} t^{\prime } = 1/(\alpha (T+B)+\beta \omega +\gamma) \end{equation} \tag{ 54 }$$

$$\begin{equation} \tilde{t}^{\prime } = 1/(\alpha U+\beta \Omega +\gamma), \end{equation} \tag{ 55 }$$

which are equivalent along the correct solution, i.e. $t^{\prime }=\tilde{t}^{\prime }$. The evolution of ω, with ω(0) = Ω(0), is obtained by

$$\begin{equation} \dot{\omega }=\sum _k{\frac{\partial \Omega }{\partial \fat r_k}}\cdot \fat v_k. \end{equation} \tag{ 56 }$$

In the presence of external perturbations the binding energy evolves according to

$$\begin{equation} \dot{B}=-\sum _k m_k {\fat v_k}\cdot \fat f_k. \end{equation} \tag{ 57 }$$

The leapfrog for the chain vectors $\fat X_k$ and $\fat V_k$ can be written as the two mappings:

$$\begin{eqnarray} &&{\fat X}(s):\hspace{19pc}\nonumber\\ &&\qquad\delta t=s/(\alpha (T+B)+\beta \omega +\gamma)\hspace{5pc} \end{eqnarray} \tag{ 58 }$$

$$\begin{equation} t=t+\delta t \end{equation} \tag{ 59 }$$

$$\begin{equation} \fat X_k\rightarrow \fat X_k+\delta t \fat V_k \end{equation} \tag{ 60 }$$

$$\begin{equation} \hspace{-19pc}{\fat V}(s): \end{equation} \tag{ 61 }$$

$$\begin{equation} \widetilde{\delta t}=s/(\alpha U+\beta \Omega +\gamma) \end{equation} \tag{ 62 }$$

$$\begin{equation} \fat V_k\rightarrow \fat V_k +\widetilde{\delta t} (\fat F_{k+1}-\fat F_k+\fat f_{k+1}-\fat f_k) \end{equation} \tag{ 63 }$$

$$\begin{equation} B\rightarrow B+\widetilde{\delta t}\sum _k(-m_k\langle \fat v_k\rangle \cdot \fat f_k) \end{equation} \tag{ 64 }$$

$$\begin{equation} \omega \rightarrow \omega +\widetilde{\delta t}\sum _k{\frac{\partial \Omega }{\partial \fat r_k}}\cdot \langle \fat v_k\rangle, \end{equation} \tag{ 65 }$$

where $\langle \fat v_k\rangle$ is the average of the initial and final $\fat v$ (obtained from the $\fat V$ according to the equations in Section 4.2.2. It is necessary to evaluate the individual velocities $\fat v_k$ because the expressions for B' and ω' in terms of the chain vector velocities $\fat V_k$ are rather cumbersome.

A leapfrog step can be written as

$$$ \fat X(h/2)\fat V(h) \fat X(h/2) $$$

and a sequence of n steps as

$$$ \fat X(h/2)\big [\Pi _{\nu =1}^{n-1} (\fat V(h) \fat X(h))\big ] \fat V(h)\fat X(h/2). $$$

This is useful with the extrapolation method when advancing the system over a total time interval of length = nh.

4.3. Time Transformation Alternatives

If one takes

$$\begin{equation} \Omega _{j}=m_im_j, \end{equation} \tag{ 66 }$$

then α = 0, β = 1, γ = 0 is mathematically equivalent to α = 1, β = γ = 0 as was shown in Mikkola & Aarseth (2002). However, numerically these are not equivalent, mainly due to roundoff errors in updating the value of Ω, and the LogH alternative is numerically more stable. However, for proper treatment of small bodies some function Ω is to be used.

For increased numerical stability in the motions of the large bodies, and smoothing of the encounters of small bodies, the recipe is α = 1, β ≠ 0 (but small) and

$$\begin{equation} \Omega _{ij}=\left\lbrace\begin{array}{@{}l@{\qquad}l@{}} {=}{\widetilde{m}}^2;& {\rm if }\quad m_i*m_j <\epsilon \: {\widetilde{m}}^2 \\ {=}0;& {\rm otherwise},\end{array}\right. \end{equation} \tag{ 67 }$$

where ${\widetilde{m}}^2=\sum _{i<j}m_im_j/(N(N-1)/2)$ is the mean mass product and ∼ 10⁻³. It may be advisable to integrate Equation (56) for ω even if β = 0, in order to force the integrator (extrapolation method!) to use short steps if $\dot{\omega }$ is large, thus giving higher precision when required. In fact, numerical experiments suggest that in most practical cases the parameters (α, β, γ) = (1, 0, 0) give the best results. Exceptions are cases with extremely large mass ratios such that the contribution of the small masses in the potential is negligible (e.g., zero masses included). This point, however, needs further investigation. One potential problem in the integration of the quantity ω is that the increments of it can be arbitrarily large if a collision of point masses occurs. In this case the roundoff errors in the value of ω are significant, but do not affect much the results if the value of β is small. One should also realize that the parameters (α, β, γ) can be changed during the integration (but not during a Bulirsch–Stoer integration step).

4.4. Velocity-Dependent Perturbations

In this section, we demonstrate the performance of the new algorithm and compare it with the celebrated KS-transformed CHAIN code of Mikkola & Aarseth (1993). We first considered a test problem consisting of one massive particle (the "black hole") with m = 1 and seven additional particles ("stars") with masses 10⁻³ ⩽ m ⩽ 10⁻⁹. Relativistic terms were not included. Initial velocities were all zero so the stars moved initially on nearly rectilinear orbits toward the black hole, but their orbits become eccentric ellipses as they experience perturbations from the other stars. We set (α, β, γ) = (1, 0, 0).

We integrated the initial values given in Table 1 including 2, 3, ..., 8 bodies and carried out integrations up to time 1000 (G = 1) using both the old KS-CHAIN code and the new AR-CHAIN. Two sets of experiments were conducted, one in which the Bulirsch–Stoer integrator automatically chose the stepsizes (in a supposedly optimal way) and the second with (iteration to) exact output times of interval Δt = 0.1.

The results are concisely summarized in Tables 2 and 3. There we give for both methods the execution times and average values of the errors. The errors are defined in terms of the Hamiltonians associated with the methods: for AR-CHAIN the average was taken over |log((T − E)/U)| ≈ δE/U and for KS-CHAIN 1.5|(T − U − E)/L| ≈ 1.5 δE/L. Here L = T + U = the Lagrangian and the factor 1.5 is included because in virial equilibrium L = 1.5U. On the other hand the system is chaotic and the solutions are, after some time, often different (as with any method) and so the results must be considered as giving just a general view of the performance.

One sees that the accuracies are comparable, with some advantage in favor of the new method. In particular, we note that the new method is faster for the case of N = 2, i.e. a binary with eccentricity e = 1, while the precision is essentially the same. Thus one can conclude that this method works equally well, and perhaps better than the celebrated KS transformation. This is a consequence of the fact that for the two-body system the LogH + leapfrog produces an exact trajectory with only a small error in the time, even for the collision orbit with eccentricity e = 1.

The accuracy of the angular momentum is typically similar to that of the energy. This is largely due to the fact that the basic algorithm (leapfrog) used in the method conserves angular momentum exactly.

However, the execution times for cases in which smaller and smaller bodies are included differ considerably in favor of the new AR-CHAIN code. This is not very surprising, since zero masses are a singularity for the KS-CHAIN, but the AR-CHAIN can also handle zero masses.

Figure 1 compares the energy conservation for the two different methods in the eight-particle integration from the initial conditions in Table 1. The system is highly chaotic and the figure shows that in this very difficult case the new method is more accurate and much faster.

Zoom In Zoom Out Reset image size

Figure 2 shows a second numerical experiment that included all PN terms up to order PN2.5. We considered the problem of periastron shift of a single star orbiting around the SMBH at the center of the Milky Way. The particle masses were m₁ = 3.5 × 10⁶ M_☉ and m₂ = 10 M_☉ = 2.85 × 10⁻⁶m₁, and the Keplerian orbit had initial semi-major axis a = 0.01 pc, similar to that of the "S" stars (Ghez et al. 2005; Eisenhauer et al. 2005). In units where G = m₁ = 1, and adopting 1 mpc = 10⁻³ pc as the length unit, the speed of light is 77.19. Three different eccentricities were tried: (0.9, 0.98, 0.99). The same values of (α, β, γ) were used as in the first experiment. The integrations were continued for 10⁶ yr or ∼20,000 orbital periods. The figure plots the orientation of the Laplace–Runge–Lenz vector every 1000 orbits; the solid (red) lines show the expected periastron advance based on the PN2.5 equations, which predict a shift each having a period of

$$\begin{equation} \Delta \phi = {6\pi GM_1\over a(1-e^2)c^2} + {3(18+e^2)\pi G^2 M_1^2\over 2a^2(1-e^2)^2c^4}, \end{equation} \tag{ 72 }$$

or (0.095, 0.45, 0.91) degrees for e = (0.9, 0.98, 0.99).

Zoom In Zoom Out Reset image size

As a third experiment, we investigated the effect of an intermediate-mass black hole (IMBH) on the motion of a star orbiting around the Milky Way SMBH. The IMBH was given a mass of 3500 M_☉ = 10⁻³M_SMBH, semi-major axis 0.1 mpc, and eccentricity 0.9. The large eccentricity greatly reduced the gravitational-wave inspiral time, by a factor ∼10³ compared with a circular orbit; inspiral required ∼3670 yr or ∼10⁵ initial orbital periods. To this binary system was added a third component of mass 10 M_☉ on an orbit having a = 8 mpc ≈ 1600 AU and e = 0.974 with respect to the center of mass of the SMBH/IMBH binary. These orbital elements are similar to those inferred for the star S0-16 (Ghez et al. 2005). The star and IMBH were coplanar with aligned angular momenta. The star completed ∼10² orbits during the time of IMBH inspiral. Figure 3 shows the evolution of the star's semi-major axis and its orbital orientation during this time, and the configuration-space trajectory is plotted in Figure 4. Initially, the apoastron distance of the IMBH is ∼0.19 mpc, similar to the periastron distance of the star, ∼0.21 mpc, so that close interactions are allowed. The primary influence was found to be on the star's semi-major axis; the eccentricity of the star's orbit changed only slightly. The star's periastron advancement was found to be well predicted by Equation (72) if the time dependence of a and e were taken into account (Figure 3). A second integration without the PN terms confirmed that the IMBH itself contributes only slightly to the star's precession rate for this configuration.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Our final numerical experiment was a one-million-year integration of a seven-body system consisting of the MW SMBH, an IMBH of mass 10⁻³M_SMBH, and the five, shortest-period S stars: S0-1, S0-2, S0-16, S0-19, and S0-20 (Ghez et al. 2005). Stellar masses were set to 15 M_☉ and the initial positions and velocities were determined at year 2000 AD using the Keplerian orbital elements given in Table 3 of Ghez et al. (2005). The IMBH orbit was assigned an initial eccentricity of 0.9 and semi-major axis of 1 mpc, compared with 4 mpc ≲ a ≲ 25 mpc for the stars. PN terms were included, causing the orbit of the IMBH to precess rapidly and to fill the annulus 0.1 mpc ≲ r ≲ 1.9 mpc; for this choice of (a, e, M_IMBH) the gravitational-wave inspiral time is ∼10⁸ yr, much greater than the length of the integration. Three of the included stars, S0-2, S0-16, and S0-19, have periastron distances that intersect the IMBH's orbital annulus and so each of these stars was able to interact closely with the IMBH, although many orbital periods were required before close encounters occurred. S0-19 was the first to achieve positive energy, at t ≈ 147, 500 AD; before ejection the star's orbit evolved toward small a and e. S0-16 was the next to escape, at t ≈ 254, 500 AD; this star moved into a highly eccentric orbit before being ejected. S0-2 remained bound to the SMBH/IMBH binary but its semi-major axis increased gradually to ∼1 pc, roughly equal to the radius of influence of the SMBH (if the latter were embedded in the Galactic bulge), so in a practical sense it too may be considered to have escaped. Figure 5 shows the evolution of a and e for the three stars. Experiments like these could be used to constrain the mass and orbital parameters of a putative IMBH near the Galactic center.

Zoom In Zoom Out Reset image size

We have demonstrated that the LogH method (with leapfrog and Bulirsch–Stoer extrapolation) is at least equally good as, and for some systems better than, the KS-CHAIN method. The optional TTL (which we use only to provide some regularization for the occasional close encounters of very small bodies) seems to have some problems that may have to do with the roundoff error when the function ω first increases then decreases by a large amount in the case of very close encounters. However, we point out that this roundoff problem is reduced almost to non-existence by including in the TTL function Ω only the interactions between small bodies and using only a small factor for this quantity in the time-transformation function. Since the TTL function derivative is explicitly integrated it automatically reduces the stepsize in the case of close approach of small bodies thus providing more careful integration of such events even if the coefficient β is negligible. This is due to the properties (sensitivity) of the Bulirsch–Stoer extrapolation method that sees any tiny irregularity in the data and modifies the stepsize in response. This property of the BS extrapolation method also helps if the leapfrog algorithm happens to evaluate the derivatives too close to the collision of some pair of bodies: the step is rejected and re-evaluated with a different stepsize. This procedure normally avoids the repetition of such a numerical accident. Finally, our code gives the option of not using time transformation at all. Even in this alternative the present code is much better than the straightforward use of center-of-mass coordinates. The reason is that the chain coordinates reduce the roundoff significantly. This is one of the great advantages of the chain structure.

D.M. was supported by grants AST-0420920 and AST-0437519 from the NSF and grant NNX07AH15G from NASA.