AN EFFICIENT CONSERVATIVE INTEGRATOR WITH A CHAIN REGULARIZATION FOR THE FEW-BODY PROBLEM

The small distance between two bodies experiencing a close encounter is represented as a difference between large numbers in straightforward formulations of the N-body problem. Thus, rounding easily becomes a significant source of error. To avoid this, we use the chain concept introduced by Mikkola & Aarseth (1993) for regularization algorithms.

In this chain method, a chain of interparticle vectors is constructed so that all the particles are included in this chain. Note that small distances are part of the chain. We begin by searching for the shortest distance, which is taken as the first piece of the chain. Next, we find the particle closest to one or the other end of the presently known part of the chain. Then, we add this particle to the end of the chain that is closer. This process is repeated until all the particles are involved. After every integration step, we check whether any non-chained vector is shorter than the smallest of the chained vectors that are in contact with one or the other end of the vector under consideration, namely, whether any triangle formed by two consecutive chained vectors has the shortest side non-chained. If this is the case, a new chain is formed. Hereafter, suppose the masses are relabeled 1, 2, ⋯, N along the chain.

We assume that ${{\boldsymbol{q}}}_{i}^{\prime }\equiv ({q}_{i[1]}^{\prime },{q}_{i[2]}^{\prime },{q}_{i[3]}^{\prime })$ is the position vector of a point with mass m_i in the barycentric frame. We also define ${{\boldsymbol{p}}}_{i}^{\prime }\equiv (p{^{\prime} }_{i[1]},p{^{\prime} }_{i[2]},p{^{\prime} }_{i[3]})$ as a momentum conjugate to ${{\boldsymbol{q}}}_{i}^{\prime }$. We set the gravitational constant equal to one for simplicity. In addition, N position vectors ${\boldsymbol{q}}{^{\prime} }_{1},\cdots ,{\boldsymbol{q}}{^{\prime} }_{N}$ and N momentum vectors ${\boldsymbol{p}}{^{\prime} }_{1},\cdots ,{\boldsymbol{p}}{^{\prime} }_{N}$ satisfy the following constraints:

$$\begin{eqnarray*}&&\displaystyle \sum _{i=1}^{N}{m}_{i}{{\boldsymbol{q}}}_{i}^{\prime }=0,\quad \displaystyle \sum _{i=1}^{N}{{\boldsymbol{p}}}_{i}^{\prime }=0.\end{eqnarray*}$$

The equations of motion in the barycentric frame are given by the Hamiltonian:

$$\begin{eqnarray}&&H=\displaystyle \frac{1}{2}\displaystyle \sum _{i=1}^{N}\displaystyle \frac{{| {{\boldsymbol{p}}}_{i}^{\prime }| }^{2}}{{m}_{i}}-\displaystyle \sum _{i=1}^{N-1}\displaystyle \sum _{j=i+1}^{N}\displaystyle \frac{{m}_{i}{m}_{j}}{| {{\boldsymbol{q}}}_{i}^{\prime }-{{\boldsymbol{q}}}_{j}^{\prime }| }.\end{eqnarray} \tag{ 1 }$$

The dynamical system corresponding to this Hamiltonian is

$$\begin{eqnarray}\displaystyle \frac{d}{{dt}}{{\boldsymbol{q}}}_{i}^{\prime } & = & \displaystyle \frac{1}{{m}_{i}}{{\boldsymbol{p}}}_{i}^{\prime },\quad \displaystyle \frac{d}{{dt}}{{\boldsymbol{p}}}_{i}^{\prime }\\ & = & {m}_{i}(\displaystyle \sum _{k=1}^{i-1}\displaystyle \frac{{m}_{k}({{\boldsymbol{q}}}_{k}^{\prime }-{{\boldsymbol{q}}}_{i}^{\prime })}{{| {{\boldsymbol{q}}}_{k}^{\prime }-{{\boldsymbol{q}}}_{i}^{\prime }| }^{3}}-\displaystyle \sum _{k=i+1}^{N}\displaystyle \frac{{m}_{k}({{\boldsymbol{q}}}_{i}^{\prime }-{{\boldsymbol{q}}}_{k}^{\prime })}{{| {{\boldsymbol{q}}}_{i}^{\prime }-{{\boldsymbol{q}}}_{k}^{\prime }| }^{3}}),\\ & & 1\leqslant i\leqslant N.\end{eqnarray} \tag{ 2 }$$

However, for two-body close encounters, we need to simultaneously use two position vectors in the barycentric frame. Therefore, the barycentric frame is not useful for computing close encounters between two masses.

The GNBP in the relative frame is much more symmetric than that in the barycentric frame. It also has a significant advantage in investigating such properties as periodic orbits and close encounters (e.g., Broucke 1975; Broucke & Boggs 1975; Mikkola & Tanikawa 1999a). It can be integrated numerically without catastrophic errors after the Levi-Civita or KS transformation (Levi-Civita 1920; Kustaanheimo & Stiefel 1965; Stiefel & Scheifele 1971).

In Section 2.2.1, we rewrite the GNBP in the relative frame. The resulting problem involves many gravitational force terms for a large number of masses N. Thus, we deform this system to reduce the number of force terms. In Section 2.2.2, we rewrite the system using the KS variables (Kustaanheimo & Stiefel 1965; Stiefel & Scheifele 1971).

2.2.1. General N-body Problem in the Relative Frame

We introduce a relative frame to consider two-body close approaches easily. We use the relative position vectors ${{\boldsymbol{q}}}_{{ij}}\equiv ({q}_{{ij}[1]},{q}_{{ij}[2]},{q}_{{ij}[3]})$ defined as

$$\begin{eqnarray}&&{{\boldsymbol{q}}}_{{ij}}={{\boldsymbol{q}}}_{i}^{\prime }-{{\boldsymbol{q}}}_{j}^{\prime },\quad 1\leqslant i\lt j\leqslant N,\end{eqnarray} \tag{ 3 }$$

and the momentum ${{\boldsymbol{p}}}_{{ij}}\equiv ({p}_{{ij}[1]},{p}_{{ij}[2]},{p}_{{ij}[3]})$ conjugate to ${{\boldsymbol{q}}}_{{ij}}$ as

$$\begin{eqnarray}&&{{\boldsymbol{p}}}_{{ij}}=\displaystyle \frac{{m}_{j}{{\boldsymbol{p}}}_{i}^{\prime }-{m}_{i}{{\boldsymbol{p}}}_{j}^{\prime }}{m},\quad 1\leqslant i\lt j\leqslant N,\end{eqnarray} \tag{ 4 }$$

where m is the total mass, ${\displaystyle \sum }_{i=1}^{N}{m}_{i}$, of the GNBP. These position and momentum vectors also satisfy the following constraints:

$$\begin{eqnarray}{{\boldsymbol{\phi }}}_{{\bf{1}}{\boldsymbol{jk}}}({\boldsymbol{q}}) & \equiv & ({\phi }_{1{jk}[1]}({\boldsymbol{q}}),{\phi }_{1{jk}[2]}({\boldsymbol{q}}),{\phi }_{1{jk}[3]}({\boldsymbol{q}}))\\ & \equiv & {{\boldsymbol{q}}}_{1j}+{{\boldsymbol{q}}}_{{jk}}-{{\boldsymbol{q}}}_{1k}=0,2\leqslant j\lt k\leqslant N,\end{eqnarray} \tag{ 5a }$$

$$\begin{eqnarray}{{\boldsymbol{\psi }}}_{{ijk}}({\boldsymbol{p}}) & \equiv & ({\psi }_{{ijk}[1]}({\boldsymbol{p}}),{\psi }_{{ijk}[2]}({\boldsymbol{p}}),{\psi }_{{ijk}[3]}({\boldsymbol{p}}))\\ & \equiv & \displaystyle \frac{{{\boldsymbol{p}}}_{{ij}}}{{m}_{i}{m}_{j}}+\displaystyle \frac{{{\boldsymbol{p}}}_{{jk}}}{{m}_{j}{m}_{k}}-\displaystyle \frac{{{\boldsymbol{p}}}_{{ik}}}{{m}_{i}{m}_{k}}=0,\\ & & 1\leqslant i\lt j\lt k\leqslant N.\end{eqnarray} \tag{ 5b }$$

Here,

$$\begin{eqnarray*}{\boldsymbol{q}} & \equiv & (\begin{array}{ccc}{q}_{\mathrm{1,2}[1]},{q}_{\mathrm{1,2}[2]},{q}_{\mathrm{1,2}[3]} & \cdots & {q}_{1,N[1]},{q}_{1,N[2]},\end{array}\\ & & {q}_{1,N[2]}| {q}_{\mathrm{2,3}[1]},{q}_{\mathrm{2,3}[2]},{q}_{\mathrm{2,3}[3]}\\ & & \begin{array}{ccc} & \cdots & {q}_{2,N[1]},{q}_{2,N[2]},{q}_{2,N[3]}| \cdots | {q}_{N-1,N[1]}\end{array},\\ & & {q}_{N-1,N[2]},{q}_{N-1,N[3]})\in {{\mathbb{R}}}^{3N(N-1)/2},\\ {\boldsymbol{p}} & \equiv & (\begin{array}{ccc}{p}_{\mathrm{1,2}[1]},{p}_{\mathrm{1,2}[2]},{p}_{\mathrm{1,2}[3]} & \cdots & {p}_{1,N[1]},{p}_{1,N[2]},\end{array}\\ & & {p}_{1,N[2]}| {p}_{\mathrm{2,3}[1]},{p}_{\mathrm{2,3}[2]},{p}_{\mathrm{2,3}[3]}\\ & & \begin{array}{ccc} & \cdots & {p}_{2,N[1]},{p}_{2,N[2]},{p}_{2,N[3]}| \cdots | {p}_{N-1,N[1]},\end{array}\\ & & {p}_{N-1,N[2]},{p}_{N-1,N[3]})\in {{\mathbb{R}}}^{3N(N-1)/2}.\end{eqnarray*}$$

To obtain the inverse transformations of Equations (3) and (4), we have to solve system (3) for the position vectors ${{\boldsymbol{q}}}_{i}^{\prime }$ and system (4) for the momentum vectors ${{\boldsymbol{p}}}_{i}^{\prime }$. Unfortunately, no vector is uniquely determined. However, if we choose

$$\begin{eqnarray}{{\boldsymbol{q}}}_{i}^{\prime } & = & \displaystyle \frac{1}{m}(\displaystyle \sum _{j=i+1}^{N}{m}_{j}{{\boldsymbol{q}}}_{{ij}}-\displaystyle \sum _{j=1}^{i-1}{m}_{j}{{\boldsymbol{q}}}_{{ji}}),\\ {{\boldsymbol{p}}}_{i}^{\prime } & = & \displaystyle \sum _{j=i+1}^{N}{{\boldsymbol{p}}}_{{ij}}-\displaystyle \sum _{j=1}^{i-1}{{\boldsymbol{p}}}_{{ji}},\ \quad 1\leqslant i\leqslant N,\end{eqnarray} \tag{ 6 }$$

these relations follow¹ Equations (3) and (4). Accordingly, we adopt Equation (6) as a transformation from the relative frame to the barycentric frame.

Substitution of transformations (3) and (4) and their time differentiations into Equation (2) yields the following system:

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dt}}{{\boldsymbol{q}}}_{{ij}}=\displaystyle \frac{m}{{m}_{i}{m}_{j}}{{\boldsymbol{p}}}_{{ij}},\quad 1\leqslant i\lt j\leqslant N,\end{eqnarray} \tag{ 7a }$$

$$\begin{eqnarray}\displaystyle \frac{d}{{dt}}{{\boldsymbol{p}}}_{1j} & = & -{m}_{1}{m}_{j}\displaystyle \frac{{{\boldsymbol{q}}}_{1j}}{{| {{\boldsymbol{q}}}_{1j}| }^{3}}+\displaystyle \sum _{k=2}^{j-1}{\boldsymbol{\mathcal{F}}}_{{kj}}({\boldsymbol{q}})\\ & & -\displaystyle \sum _{k=j+1}^{N}{\boldsymbol{\mathcal{F}}}_{{jk}}({\boldsymbol{q}}),\quad 2\leqslant j\leqslant N,\end{eqnarray} \tag{ 7b }$$

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dt}}{{\boldsymbol{p}}}_{{ij}}=-{m}_{i}{m}_{j}\displaystyle \frac{{{\boldsymbol{q}}}_{{ij}}}{{| {{\boldsymbol{q}}}_{{ij}}| }^{3}}-{\boldsymbol{\mathcal{F}}}_{{ij}}({\boldsymbol{q}}),\quad 2\leqslant i\lt j\leqslant N,\end{eqnarray} \tag{ 7c }$$

where ${\boldsymbol{\mathcal{F}}}_{{ij}}({\boldsymbol{q}})$ is a function of a vector ${\boldsymbol{q}}$ defined by

$$\begin{eqnarray*}{\boldsymbol{\mathcal{F}}}_{{ij}}({\boldsymbol{q}}) & \equiv & -\displaystyle \sum _{k=1}^{i-1}\displaystyle \frac{{m}_{i}{m}_{j}{m}_{k}}{m}(\displaystyle \frac{{{\boldsymbol{q}}}_{{ki}}}{{| {{\boldsymbol{q}}}_{{ki}}| }^{3}}+\displaystyle \frac{{{\boldsymbol{q}}}_{{ij}}}{{| {{\boldsymbol{q}}}_{{ij}}| }^{3}}-\displaystyle \frac{{{\boldsymbol{q}}}_{{kj}}}{{| {{\boldsymbol{q}}}_{{kj}}| }^{3}})\\ & & +\displaystyle \sum _{k=i+1}^{j-1}\displaystyle \frac{{m}_{i}{m}_{j}{m}_{k}}{m}(\displaystyle \frac{{{\boldsymbol{q}}}_{{ik}}}{{| {{\boldsymbol{q}}}_{{ik}}| }^{3}}+\displaystyle \frac{{{\boldsymbol{q}}}_{{kj}}}{{| {{\boldsymbol{q}}}_{{kj}}| }^{3}}-\displaystyle \frac{{{\boldsymbol{q}}}_{{ij}}}{{| {{\boldsymbol{q}}}_{{ij}}| }^{3}})\\ & & -\displaystyle \sum _{k=j+1}^{N}\displaystyle \frac{{m}_{i}{m}_{j}{m}_{k}}{m}(\displaystyle \frac{{{\boldsymbol{q}}}_{{ij}}}{{| {{\boldsymbol{q}}}_{{ij}}| }^{3}}+\displaystyle \frac{{{\boldsymbol{q}}}_{{jk}}}{{| {{\boldsymbol{q}}}_{{jk}}| }^{3}}-\displaystyle \frac{{{\boldsymbol{q}}}_{{ik}}}{{| {{\boldsymbol{q}}}_{{ik}}| }^{3}}),\\ & & \quad 2\leqslant i\lt j\leqslant N.\end{eqnarray*}$$

We can regard the system composed of Equations (5(a)) and (7) as the GNBP in the relative frame.

In combination with Equation (6), the Hamiltonian (1) leads to

$$\begin{eqnarray*}H & = & \displaystyle \sum _{i=1}^{N-1}\displaystyle \sum _{j=i+1}^{N}(\displaystyle \frac{m}{2{m}_{i}{m}_{j}}{| {{\boldsymbol{p}}}_{{ij}}| }^{2}-\displaystyle \frac{{m}_{i}{m}_{j}}{| {{\boldsymbol{q}}}_{{ij}}| })\\ & & +\displaystyle \frac{1}{2}\displaystyle \sum _{i=1}^{N-2}\displaystyle \sum _{j=i+1}^{N-1}\;\displaystyle \sum _{k=j+1}^{N}{m}_{i}{m}_{j}{m}_{k}{| {{\boldsymbol{\psi }}}_{{ijk}}({\boldsymbol{p}})| }^{2}.\end{eqnarray*}$$

However, it does not yield system (7), so it cannot be regarded as a Hamiltonian. System (7) is governed by the following Hamiltonian:

$$\begin{eqnarray}{H}_{\mathrm{rel}} & \equiv & \displaystyle \sum _{i=1}^{N-1}\displaystyle \sum _{j=i+1}^{N}(\displaystyle \frac{m}{2{m}_{i}{m}_{j}}{| {{\boldsymbol{p}}}_{{ij}}| }^{2}-\displaystyle \frac{{m}_{i}{m}_{j}}{| {{\boldsymbol{q}}}_{{ij}}| })\\ & & +\displaystyle \sum _{j=2}^{N-1}\displaystyle \sum _{k=j+1}^{N}{{\boldsymbol{\phi }}}_{1{jk}}({\boldsymbol{q}})\cdot {{\boldsymbol{\lambda }}}_{{jk}},\end{eqnarray} \tag{ 8 }$$

where ${{\boldsymbol{\lambda }}}_{{jk}}=({\lambda }_{{jk}[1]},{\lambda }_{{jk}[2]},{\lambda }_{{jk}[3]})$ are Lagrange multipliers.² The values of H and ${H}_{\mathrm{rel}}$ coincide with

$$\begin{eqnarray}&&{h}_{\mathrm{rel}}\equiv \displaystyle \sum _{i=1}^{N-1}\displaystyle \sum _{j=i+1}^{N}(\displaystyle \frac{m}{2{m}_{i}{m}_{j}}{| {{\boldsymbol{p}}}_{{ij}}| }^{2}-\displaystyle \frac{{m}_{i}{m}_{j}}{| {{\boldsymbol{q}}}_{{ij}}| }).\end{eqnarray} \tag{ 9 }$$

2.2.2. General N-body Problem with KS Variables

In the numerical integration of the GNBP, multi-body close encounters result in serious numerical difficulties due to the errors associated with singularities in the GNBP. KS regularization is a standard technique for removing the singularities. It combines a time regularization with the KS transformation (Kustaanheimo & Stiefel 1965). In this subsection, using the KS variables (Kustaanheimo & Stiefel 1965), we rewrite the GNBP in Equation (7).

We apply the KS transformation (Kustaanheimo & Stiefel 1965) to the vectors ${\bar{{\boldsymbol{q}}}}_{{ij}}\equiv ({{\boldsymbol{q}}}_{{ij}}\ | \ 0)=({q}_{{ij}[1]},{q}_{{ij}[2]},{q}_{{ij}[3]},0)$, obtaining the vectors ${{\boldsymbol{Q}}}_{{ij}}\equiv ({Q}_{{ij}[1]},{Q}_{{ij}[2]},{Q}_{{ij}[3]},{Q}_{{ij}[4]})$. The relations between ${\bar{{\boldsymbol{q}}}}_{{ij}}$ and ${{\boldsymbol{Q}}}_{{ij}}$ are given by

$$\begin{eqnarray}&&{\bar{{\boldsymbol{q}}}}_{{ij}}={{\boldsymbol{Q}}}_{{ij}}{\boldsymbol{K}}{({{\boldsymbol{Q}}}_{{ij}})}^{\top },\quad 1\leqslant i\lt j\leqslant N,\end{eqnarray} \tag{ 10 }$$

where the KS matrix (Kustaanheimo & Stiefel 1965) is defined as

$$\begin{eqnarray*}{\boldsymbol{K}}({{\boldsymbol{Q}}}_{{ij}})\equiv [\begin{array}{cccc}{Q}_{{ij}[1]} & -{Q}_{{ij}[2]} & -{Q}_{{ij}[3]} & {Q}_{{ij}[4]}\\ {Q}_{{ij}[2]} & {Q}_{{ij}[1]} & -{Q}_{{ij}[4]} & -{Q}_{{ij}[3]}\\ {Q}_{{ij}[3]} & {Q}_{{ij}[4]} & {Q}_{{ij}[1]} & {Q}_{{ij}[2]}\\ {Q}_{{ij}[4]} & -{Q}_{{ij}[3]} & {Q}_{{ij}[2]} & -{Q}_{{ij}[1]}\end{array}].\end{eqnarray*}$$

New momentum vectors ${{\boldsymbol{P}}}_{{ij}}\equiv ({P}_{{ij}[1]},{P}_{{ij}[2]},{P}_{{ij}[3]},{P}_{{ij}[4]})$ are also related to ${\bar{{\boldsymbol{p}}}}_{{ij}}\equiv ({{\boldsymbol{p}}}_{{ij}}\ | \ 0)=({p}_{{ij}[1]},{p}_{{ij}[2]},{p}_{{ij}[3]},0)$ by the equations

$$\begin{eqnarray}&&{\bar{{\boldsymbol{p}}}}_{{ij}}=\displaystyle \frac{1}{2{| {{\boldsymbol{Q}}}_{{ij}}| }^{2}}{{\boldsymbol{P}}}_{{ij}}{\boldsymbol{K}}{({{\boldsymbol{Q}}}_{{ij}})}^{\top },1\leqslant i\lt j\leqslant N,\end{eqnarray} \tag{ 11 }$$

where

$$\begin{eqnarray*} & & {P}_{{ij}[1]}{Q}_{{ij}[4]}-{P}_{{ij}[2]}{Q}_{{ij}[3]}+{P}_{{ij}[3]}{Q}_{{ij}[2]}\\ & & -{P}_{{ij}[4]}{Q}_{{ij}[1]}=0,\quad 1\leqslant i\lt j\leqslant N.\end{eqnarray*}$$

Theoretically, we can use the solution of relations (10) and (11) to obtain ${{\boldsymbol{Q}}}_{{ij}}$ and ${{\boldsymbol{P}}}_{{ij}}$ from ${{\boldsymbol{q}}}_{{ij}}$ and ${{\boldsymbol{p}}}_{{ij}}$, respectively.

In numerical computation, to avoid cancellation of significant digits, we compute ${{\boldsymbol{Q}}}_{{ij}}$ for $1\leqslant i\lt j\leqslant N$ as follows:

$$\begin{eqnarray} & & {{\boldsymbol{Q}}}_{{ij}}\equiv \{\begin{array}{cc}({u}_{{ij}[1]},{u}_{{ij}[2]},{u}_{{ij}[3]},{u}_{{ij}[4]}) & {\rm{if}}\ {q}_{{ij}[1]}\geqslant 0,\\ ({u}_{{ij}[2]},{u}_{{ij}[1]},{u}_{{ij}[4]},{u}_{{ij}[3]}) & {\rm{if}}\ {q}_{{ij}[1]}\lt 0,\end{array}\end{eqnarray} \tag{ 12 }$$

where

$$\begin{eqnarray*}{u}_{{ij}[1]} & = & {(\displaystyle \frac{| {{\boldsymbol{q}}}_{{ij}}| +| {q}_{{ij}[1]}| }{2})}^{1/2},\\ {u}_{{ij}[2]} & = & \displaystyle \frac{{q}_{{ij}[2]}}{2{u}_{{ij}[1]}},\ {u}_{{ij}[3]}=\displaystyle \frac{{q}_{{ij}[3]}}{2{u}_{{ij}[1]}},\\ {u}_{{ij}[4]} & = & 0,1\leqslant i\lt j\leqslant N.\end{eqnarray*}$$

The KS momenta ${{\boldsymbol{P}}}_{{ij}}$ are obtained as

$$\begin{eqnarray*}&&{{\boldsymbol{P}}}_{{ij}}=2{\bar{{\boldsymbol{p}}}}_{{ij}}{\boldsymbol{K}}({{\boldsymbol{Q}}}_{{ij}}),\quad 1\leqslant i\lt j\leqslant N.\end{eqnarray*}$$

Substitution of transformations (10) and (11) and their time differentiations into Equation (7) yields

$$\begin{eqnarray}&&{{\boldsymbol{F}}}_{{ij}}=0,1\leqslant i\lt j\leqslant N,\end{eqnarray} \tag{ 13a }$$

$$\begin{eqnarray}&&{{\boldsymbol{G}}}_{1j}-\displaystyle \sum _{k=2}^{j-1}{\bar{{\boldsymbol{\lambda }}}}_{{kj}}+\displaystyle \sum _{k=j+1}^{N}{\bar{{\boldsymbol{\lambda }}}}_{{jk}}=0,2\leqslant j\leqslant N,\end{eqnarray} \tag{ 13b }$$

$$\begin{eqnarray}&&{{\boldsymbol{G}}}_{{ij}}+{\bar{{\boldsymbol{\lambda }}}}_{{ij}}=0,2\leqslant i\lt j\leqslant N,\end{eqnarray} \tag{ 13c }$$

where

$$\begin{eqnarray}&&{{\boldsymbol{F}}}_{{ij}}\equiv \displaystyle \frac{d}{{dt}}{{\boldsymbol{Q}}}_{{ij}}-\displaystyle \frac{m}{4{m}_{i}{m}_{j}}\displaystyle \frac{{{\boldsymbol{P}}}_{{ij}}}{{| {{\boldsymbol{Q}}}_{{ij}}| }^{2}},1\leqslant i\lt j\leqslant N,\end{eqnarray} \tag{ 14a }$$

$$\begin{eqnarray}{{\boldsymbol{G}}}_{{ij}} & \equiv & \displaystyle \frac{1}{2{| {{\boldsymbol{Q}}}_{{ij}}| }^{2}}(\displaystyle \frac{d}{{dt}}{{\boldsymbol{P}}}_{{ij}}-(\displaystyle \frac{m}{4{m}_{i}{m}_{j}}{| {{\boldsymbol{P}}}_{{ij}}| }^{2}-2{m}_{i}{m}_{j})\displaystyle \frac{{{\boldsymbol{Q}}}_{{ij}}}{{| {{\boldsymbol{Q}}}_{{ij}}| }^{4}})\\ & & {{\boldsymbol{K}}}^{\top }({{\boldsymbol{Q}}}_{{ij}}),1\leqslant i\lt j\leqslant N,\end{eqnarray} \tag{ 14b }$$

$$\begin{eqnarray}&&{\bar{{\boldsymbol{\lambda }}}}_{{ij}}\equiv -\displaystyle \sum _{k=1}^{i-1}\displaystyle \frac{{m}_{i}{m}_{j}{m}_{k}}{m}\\ &&\times (\displaystyle \frac{{{\boldsymbol{Q}}}_{{ki}}}{{| {{\boldsymbol{Q}}}_{{ki}}| }^{6}}{{\boldsymbol{K}}}^{\top }({{\boldsymbol{Q}}}_{{ki}})+\displaystyle \frac{{{\boldsymbol{Q}}}_{{ij}}}{{| {{\boldsymbol{Q}}}_{{ij}}| }^{6}}{{\boldsymbol{K}}}^{\top }({{\boldsymbol{Q}}}_{{ij}})-\displaystyle \frac{{{\boldsymbol{Q}}}_{{kj}}}{{| {{\boldsymbol{Q}}}_{{kj}}| }^{6}}{{\boldsymbol{K}}}^{\top }({{\boldsymbol{Q}}}_{{kj}}))\\ &&+\displaystyle \sum _{k=i+1}^{j-1}\displaystyle \frac{{m}_{i}{m}_{j}{m}_{k}}{m}\\ &&\times (\displaystyle \frac{{{\boldsymbol{Q}}}_{{ik}}}{{| {{\boldsymbol{Q}}}_{{ik}}| }^{6}}{{\boldsymbol{K}}}^{\top }({{\boldsymbol{Q}}}_{{ik}})+\displaystyle \frac{{{\boldsymbol{Q}}}_{{kj}}}{{| {{\boldsymbol{Q}}}_{{kj}}| }^{6}}{{\boldsymbol{K}}}^{\top }({{\boldsymbol{Q}}}_{{kj}})-\displaystyle \frac{{{\boldsymbol{Q}}}_{{ij}}}{{| {{\boldsymbol{Q}}}_{{ij}}| }^{6}}{{\boldsymbol{K}}}^{\top }({{\boldsymbol{Q}}}_{{ij}}))\\ &&-\displaystyle \sum _{k=j+1}^{N}\displaystyle \frac{{m}_{i}{m}_{j}{m}_{k}}{m}\\ &&\times (\displaystyle \frac{{{\boldsymbol{Q}}}_{{ij}}}{{| {{\boldsymbol{Q}}}_{{ij}}| }^{6}}{{\boldsymbol{K}}}^{\top }({{\boldsymbol{Q}}}_{{ij}})+\displaystyle \frac{{{\boldsymbol{Q}}}_{{jk}}}{{| {{\boldsymbol{Q}}}_{{jk}}| }^{6}}{{\boldsymbol{K}}}^{\top }({{\boldsymbol{Q}}}_{{jk}})-\displaystyle \frac{{{\boldsymbol{Q}}}_{{ik}}}{{| {{\boldsymbol{Q}}}_{{ik}}| }^{6}}{{\boldsymbol{K}}}^{\top }({{\boldsymbol{Q}}}_{{ik}})),\\ &&2\leqslant i\lt j\leqslant N.\end{eqnarray} \tag{ 14c }$$

In addition, in combination with Equation (10), the constraints in Equation (5(a)) lead to

$$\begin{eqnarray}{{\boldsymbol{\Phi }}}_{1{jk}}({\boldsymbol{Q}}) & \equiv & {[\begin{array}{c}{{\rm{\Phi }}}_{1{jk}[1]}({\boldsymbol{Q}})\\ {{\rm{\Phi }}}_{1{jk}[2]}({\boldsymbol{Q}})\\ {{\rm{\Phi }}}_{1{jk}[3]}({\boldsymbol{Q}})\\ 0\end{array}]}^{\top }\\ & \equiv & {{\boldsymbol{Q}}}_{1j}{\boldsymbol{K}}{({{\boldsymbol{Q}}}_{1j})}^{\top }+{{\boldsymbol{Q}}}_{{jk}}{\boldsymbol{K}}{({{\boldsymbol{Q}}}_{{jk}})}^{\top }\\ & & -{{\boldsymbol{Q}}}_{1k}{\boldsymbol{K}}{({{\boldsymbol{Q}}}_{1k})}^{\top }=0,2\leqslant j\lt k\leqslant N,\end{eqnarray} \tag{ 15 }$$

where $\begin{array}{ccccc}{\boldsymbol{Q}}\equiv ({{\boldsymbol{Q}}}_{\mathrm{1,2}}, & \cdots & {{\boldsymbol{Q}}}_{1,N},| {{\boldsymbol{Q}}}_{\mathrm{2,3}}, & \cdots & {{\boldsymbol{Q}}}_{2,N},\end{array}$ $| \cdots | {{\boldsymbol{Q}}}_{N-1,N})$ $\in {{\mathbb{R}}}^{2N(N-1)}$ is a position vector. The system composed of Equations (13) and (15) represents the motion of the GNBP, which is described by the KS variables and the Lagrange multipliers ${\bar{{\boldsymbol{\lambda }}}}_{{jk}}\equiv ({\bar{{\boldsymbol{\lambda }}}}_{{jk}}| 0)=({\lambda }_{{jk}[1]},{\lambda }_{{jk}[2]},{\lambda }_{{jk}[3]},0)\in {{\mathbb{R}}}^{4}$. This system is governed by the following Hamiltonian:

$$\begin{eqnarray}{H}_{{\rm{LC}}} & \equiv & \displaystyle \sum _{i=1}^{N-1}\displaystyle \sum _{j=i+1}^{N}(\displaystyle \frac{m}{8{m}_{i}{m}_{j}}\displaystyle \frac{{| {{\boldsymbol{P}}}_{{ij}}| }^{2}}{{| {{\boldsymbol{Q}}}_{{ij}}| }^{2}}-\displaystyle \frac{{m}_{i}{m}_{j}}{{| {{\boldsymbol{Q}}}_{{ij}}| }^{2}})\\ & & +\displaystyle \sum _{j=2}^{N-1}\displaystyle \sum _{k=j+1}^{N}{{\boldsymbol{\Phi }}}_{1{jk}}({\boldsymbol{Q}})\cdot {\bar{{\boldsymbol{\lambda }}}}_{{jk}},\end{eqnarray} \tag{ 16 }$$

which is obtained by substituting Equations (10) and (11) into ${H}_{{\rm{r}}{\rm{e}}{\rm{l}}}$ as defined by Equation (8). In addition, using Equation (11), we obtain the relation

$$\begin{eqnarray*}&&{| {{\boldsymbol{p}}}_{{ij}}| }^{2}=\displaystyle \frac{{| {{\boldsymbol{P}}}_{{ij}}| }^{2}}{4{| {{\boldsymbol{Q}}}_{{ij}}| }^{2}},\quad 1\leqslant i\lt j\leqslant N.\end{eqnarray*}$$

Using Equation (13(c)), we rewrite Equation (13(b)) as the following identities:

$$\begin{eqnarray}&&\displaystyle \sum _{i=1}^{j-1}{{\boldsymbol{G}}}_{{ij}}-\displaystyle \sum _{i=j+1}^{N}{{\boldsymbol{G}}}_{{ji}}=0,2\leqslant j\leqslant N.\end{eqnarray} \tag{ 17 }$$

The new system composed of Equations (13a), (15), and (17) describes the same motion as the system composed of Equations (13) and (15). The number of dependent variables ${P}_{{ij}[1]}$, ${P}_{{ij}[2]}$, ${P}_{{ij}[3]}$, ${P}_{{ij}[4]}$, ${Q}_{{ij}[1]}$, ${Q}_{{ij}[2]}$, ${Q}_{{ij}[3]}$, and ${Q}_{{ij}[4]}$ in the new system is $4N(N-1)$ larger than $6N$, which is the number of dependent variables $p{^{\prime} }_{i[1]}$, $p{^{\prime} }_{i[2]}$, $p{^{\prime} }_{i[3]}$, ${q}_{i[1]}^{\prime }$, ${q}_{i[2]}^{\prime }$, and ${q}_{i[3]}^{\prime }$ in system (2). Actually, as the total number of masses N increases, the number of equations in the new system increases more rapidly. Therefore, the computational cost of the new system is very high for large N. We call this system the redundant general N-body problem (RGNBP).

${H}_{{\rm{LC}}}$ is conserved by the RGNBP as well as by the system composed of Equations (13) and (15). Because of Equation (15), the value of ${H}_{{\rm{LC}}}$ defined by Equation (16) is equivalent to

$$\begin{eqnarray}&&{h}_{{\rm{LC}}}=\displaystyle \sum _{i=1}^{N-1}\displaystyle \sum _{j=i+1}^{N}(\displaystyle \frac{m}{8{m}_{i}{m}_{j}}\displaystyle \frac{{| {{\boldsymbol{P}}}_{{ij}}| }^{2}}{{| {{\boldsymbol{Q}}}_{{ij}}| }^{2}}-\displaystyle \frac{{m}_{i}{m}_{j}}{| {{\boldsymbol{Q}}}_{{ij}}| }).\end{eqnarray} \tag{ 18 }$$

Using Equations (10) and (11), the value of the Hamiltonian ${h}_{\mathrm{rel}}$ in Equation (9) reduces to ${h}_{{\rm{LC}}}$, so ${h}_{{\rm{LC}}}={h}_{\mathrm{rel}}$.

Further, ${h}_{{\rm{LC}}}$ equals the value of H as defined by Equation (1) because ${h}_{\mathrm{rel}}$ as defined by Equation (9) is transformed to ${h}_{{\rm{LC}}}$ using Equations (10) and (11), and ${h}_{\mathrm{rel}}$ is the value of H as described in Section 2.2.1. In this section, we rewrite the RGNBP using only ${{\boldsymbol{P}}}_{k,k+1}$ and ${{\boldsymbol{Q}}}_{k,k+1}$ $(1\leqslant k\leqslant N-1)$ to reduce the redundancy of the problem, which incurs a high computational cost. The vectors ${{\boldsymbol{Q}}}_{k,k+1}$ and ${{\boldsymbol{P}}}_{k,k+1}$ are related to the chained position vectors ${{\boldsymbol{q}}}_{k,k+1}$ and their momenta ${{\boldsymbol{p}}}_{k,k+1}$ conjugate to ${{\boldsymbol{q}}}_{k,k+1}$, respectively.

First, we express the vectors ${{\boldsymbol{Q}}}_{{ij}}$ related to the non-chained position vectors ${{\boldsymbol{q}}}_{{ij}}$ $(1\leqslant i\lt i+2\leqslant j\leqslant N)$ as functions of ${{\boldsymbol{Q}}}_{k,k+1}$ related to the chained position vectors ${{\boldsymbol{q}}}_{k,k+1}$ $(1\leqslant k\leqslant N-1)$. The constraints in Equation (15) restrict the possible positions of the RGNBP to the $2N(N-1)-2(N-1)(N-2)=4(N-1)$-dimensional configuration manifold

$$\begin{eqnarray}&&{\mathcal{Q}}=\{{\boldsymbol{Q}}\in {{\mathbb{R}}}^{2N(N-1)}\ | \ {\boldsymbol{\Phi }}({\boldsymbol{Q}})={0}_{1\times 2(N-1)(N-2)}\},\end{eqnarray} \tag{ 19 }$$

where $\begin{array}{ccccc}{\boldsymbol{\Phi }}({\boldsymbol{Q}})=({{\boldsymbol{\Phi }}}_{\mathrm{1,2,3}}({\boldsymbol{Q}}), & \cdots & {{\boldsymbol{\Phi }}}_{\mathrm{1,2},N}({\boldsymbol{Q}}),| {{\boldsymbol{\Phi }}}_{\mathrm{1,3,4}}({\boldsymbol{Q}}), & \cdots & {{\boldsymbol{\Phi }}}_{\mathrm{1,3},N}({\boldsymbol{Q}}),\end{array}$ $| \cdots | {{\boldsymbol{\Phi }}}_{1,N-1,N}({\boldsymbol{Q}}))\in {{\mathbb{R}}}^{2(N-1)(N-2)}$ is a constraint function vector. Accordingly, the set of $4(N-1)$ variables $\{{Q}_{\mathrm{1,2}[1]},{Q}_{\mathrm{1,2}[2]},{Q}_{\mathrm{1,2}[3]},$ $\begin{array}{cccc}{Q}_{\mathrm{1,2}[4]}, & {Q}_{\mathrm{2,3}[1]},\cdots ,{Q}_{\mathrm{2,3}[4]}, & \cdots & {Q}_{N-1,N[1]},\cdots ,{Q}_{N-1,N[4]}\end{array}\}$, which corresponds to the $(N-1)$ vectors ${{\boldsymbol{Q}}}_{k,k+1}$ $(1\leqslant k\leqslant N-1)$, acts as a basis for the manifold ${Q}$.

In addition to the constraints in Equation (5(a)), the relative position vectors ${{\boldsymbol{q}}}_{{ij}}$ satisfy

$$\begin{eqnarray}&&{{\boldsymbol{q}}}_{{ij}}=\displaystyle \sum _{k=i}^{j-1}{{\boldsymbol{q}}}_{k,k+1},\quad 1\leqslant i\lt i+2\leqslant j\leqslant N.\end{eqnarray} \tag{ 20 }$$

Substitution of Equation (10) into the right-hand side of Equation (20) yields

$$\begin{eqnarray}{{\boldsymbol{q}}}_{{ij}} & = & {\tilde{{\boldsymbol{q}}}}_{{ij}}\equiv {[\begin{array}{c}{\tilde{q}}_{{ij}[1]}\\ {\tilde{q}}_{{ij}[2]}\\ {\tilde{q}}_{{ij}[3]}\end{array}]}^{\top }\\ & = & {[\begin{array}{c}{\displaystyle \sum }_{k=i}^{j-1}({({Q}_{k,k+1[1]})}^{2}-{({Q}_{k,k+1[2]})}^{2}\\ \quad -\;{({Q}_{k,k+1[3]})}^{2}+{({Q}_{k,k+1[4]})}^{2})\\ 2{\displaystyle \sum }_{k=i}^{j-1}({Q}_{k,k+1[1]}{Q}_{k,k+1[2]}\\ \quad -\;{Q}_{k,k+1[3]}{Q}_{k,k+1[4]})\\ 2{\displaystyle \sum }_{k=i}^{j-1}({Q}_{k,k+1[1]}{Q}_{k,k+1[3]}\\ \quad +{Q}_{k,k+1[2]}{Q}_{k,k+1[4]})\end{array}]}^{\top }.\end{eqnarray} \tag{ 21 }$$

We subsequently substitute Equation (21) into (12). Then, we obtain

$$\begin{eqnarray}{{\boldsymbol{Q}}}_{{ij}} & = & {\tilde{{\boldsymbol{Q}}}}_{{ij}}\equiv ({\tilde{Q}}_{{ij}[1]},{\tilde{Q}}_{{ij}[2]},{\tilde{Q}}_{{ij}[3]},{\tilde{Q}}_{{ij}[4]})\\ & = & \{\begin{array}{cc}({\tilde{u}}_{{ij}[1]},{\tilde{u}}_{{ij}[2]},{\tilde{u}}_{{ij}[3]},{\tilde{u}}_{{ij}[4]}) & {\rm{if}}\ {\tilde{q}}_{{ij}[1]}\geqslant 0,\\ ({\tilde{u}}_{{ij}[2]},{\tilde{u}}_{{ij}[1]},{\tilde{u}}_{{ij}[4]},{\tilde{u}}_{{ij}[3]}) & {\rm{if}}\ {\tilde{q}}_{{ij}[1]}\lt 0,\end{array}\end{eqnarray} \tag{ 22 }$$

where

$$\begin{eqnarray*}{\tilde{u}}_{{ij}[1]} & = & {(\displaystyle \frac{| {\tilde{{\boldsymbol{q}}}}_{{ij}}| +| {\tilde{q}}_{{ij}[1]}| }{2})}^{1/2},\\ {\tilde{u}}_{{ij}[2]} & = & \displaystyle \frac{{\tilde{q}}_{{ij}[2]}}{2{\tilde{u}}_{{ij}[1]}},\;{\tilde{u}}_{{ij}[3]}=\displaystyle \frac{{\tilde{q}}_{{ij}[3]}}{2{\tilde{u}}_{{ij}[1]}},\\ {\tilde{u}}_{{ij}[4]} & = & 0,1\leqslant i\lt j\leqslant N.\end{eqnarray*}$$

Further, we expand Equation (20) to the extended relations, ${\bar{{\boldsymbol{q}}}}_{{ij}}={\displaystyle \sum }_{k=i}^{j-1}{\bar{{\boldsymbol{q}}}}_{k,k+1}$ $(1\leqslant i\lt i+2\leqslant j\leqslant N)$, and subsequently substitute Equation (22) into these extended relations. Then, we have

$$\begin{eqnarray}{{\boldsymbol{Q}}}_{{ij}}{\boldsymbol{K}}{({{\boldsymbol{Q}}}_{{ij}})}^{\top } & = & \displaystyle \sum _{k=i}^{j-1}{{\boldsymbol{Q}}}_{k,k+1}{\boldsymbol{K}}{({{\boldsymbol{Q}}}_{k,k+1})}^{\top },\\ & & \quad 1\leqslant i\lt i+2\leqslant j\leqslant N.\end{eqnarray} \tag{ 23 }$$

Next, we write the vectors ${{\boldsymbol{P}}}_{{ij}}$ $(1\leqslant i\lt i+2\leqslant j\leqslant N)$ related to the momentum vectors ${{\boldsymbol{p}}}_{{ij}}$ as functions of ${{\boldsymbol{Q}}}_{k,k+1}$ and ${{\boldsymbol{P}}}_{k,k+1}$ $(1\leqslant k\leqslant N-1)$ related to the momentum vectors ${{\boldsymbol{p}}}_{k,k+1}$. The momenta of the RGNBP, ${{\boldsymbol{P}}}_{{jk}}$, satisfy the following constraints:

$$\begin{eqnarray}{{\boldsymbol{\Psi }}}_{1{jk}}({\boldsymbol{P}},{\boldsymbol{Q}}) & \equiv & \displaystyle \frac{1}{{m}_{1}{m}_{j}}\displaystyle \frac{{{\boldsymbol{P}}}_{1j}{\boldsymbol{K}}{({{\boldsymbol{Q}}}_{1j})}^{\top }}{{| {{\boldsymbol{Q}}}_{1j}| }^{2}}\\ & & +\displaystyle \frac{1}{{m}_{j}{m}_{k}}\displaystyle \frac{{{\boldsymbol{P}}}_{{jk}}{\boldsymbol{K}}{({{\boldsymbol{Q}}}_{{jk}})}^{\top }}{{| {{\boldsymbol{Q}}}_{{jk}}| }^{2}}\\ & & -\displaystyle \frac{1}{{m}_{1}{m}_{k}}\displaystyle \frac{{{\boldsymbol{P}}}_{1k}{\boldsymbol{K}}{({{\boldsymbol{Q}}}_{1k})}^{\top }}{{| {{\boldsymbol{Q}}}_{1k}| }^{2}}=0,\\ & & 2\leqslant j\lt k\leqslant N,\end{eqnarray} \tag{ 24 }$$

which is obtained by substituting Equation (13(a)) into the time derivative of ${{\boldsymbol{\Phi }}}_{1{jk}}({\boldsymbol{Q}})$ in Equation (15). Here, $\begin{array}{ccccc}{\boldsymbol{P}}=({{\boldsymbol{P}}}_{\mathrm{1,2}}, & \cdots & {{\boldsymbol{P}}}_{1,N},| {{\boldsymbol{P}}}_{\mathrm{2,3}}, & \cdots & {{\boldsymbol{P}}}_{2,N},\end{array}$ $| \cdots | {{\boldsymbol{P}}}_{N-1,N})\in {{\mathbb{R}}}^{2N(N-1)}$ is a momentum vector. The constraints in Equation (24) restrict the possible momenta of the RGNBP to the $2N(N-1)-2(N-1)(N-2)=4(N-1)$-dimensional manifold

$$\begin{eqnarray*}&&{\mathcal{P}}=\{{\boldsymbol{P}}\in {{\mathbb{R}}}^{2N(N-1)}\ | \ {\boldsymbol{\Psi }}({\boldsymbol{P}},{\boldsymbol{Q}})={0}_{1\times 2(N-1)(N-2)}\},\end{eqnarray*}$$

where

$$\begin{eqnarray*}{\boldsymbol{\Psi }}({\boldsymbol{P}},{\boldsymbol{Q}}) & = & (\begin{array}{ccc}{{\boldsymbol{\Psi }}}_{\mathrm{1,2,3}}({\boldsymbol{P}},{\boldsymbol{Q}}), & \cdots & {{\boldsymbol{\Psi }}}_{\mathrm{1,2},N}({\boldsymbol{P}},{\boldsymbol{Q}}),\end{array}\\ & & | \begin{array}{ccc}{{\boldsymbol{\Psi }}}_{\mathrm{1,3,4}}({\boldsymbol{P}},{\boldsymbol{Q}}), & \cdots & {{\boldsymbol{\Psi }}}_{\mathrm{1,3},N}({\boldsymbol{P}},{\boldsymbol{Q}}),\end{array}\\ & & | \cdots | {{\boldsymbol{\Psi }}}_{1,N-1,N}({\boldsymbol{P}},{\boldsymbol{Q}}))\\ & & \in {{\mathbb{R}}}^{2(N-1)(N-2)}\end{eqnarray*}$$

is a constraint function vector. Thus, the $(N-1)$ vectors ${{\boldsymbol{P}}}_{k,k+1}$ $(1\leqslant k\leqslant N-1)$ are a basis for the manifold ${\mathcal{P}}$. Each vector ${{\boldsymbol{P}}}_{{ij}}\in {{\mathbb{R}}}^{2(N-1)(N-2)}$ $(1\leqslant i\lt i+2\leqslant j\leqslant N)$ related to a non-chained momentum vector ${{\boldsymbol{p}}}_{{ij}}$ can be uniquely expressed as follows:

$$\begin{eqnarray}\displaystyle \frac{{{\boldsymbol{P}}}_{{ij}}{\boldsymbol{K}}{({{\boldsymbol{Q}}}_{{ij}})}^{\top }}{{| {{\boldsymbol{Q}}}_{{ij}}| }^{2}} & = & {m}_{i}{m}_{j}\\ & & \times \displaystyle \sum _{k=i}^{j-1}\displaystyle \frac{1}{{m}_{k}{m}_{k+1}}\displaystyle \frac{{{\boldsymbol{P}}}_{k,k+1}{\boldsymbol{K}}{({{\boldsymbol{Q}}}_{k,k+1})}^{\top }}{{| {{\boldsymbol{Q}}}_{k,k+1}| }^{2}},\\ & & 1\leqslant i\lt i+2\leqslant j\leqslant N.\end{eqnarray} \tag{ 25 }$$

We solve Equation (25) for ${{\boldsymbol{P}}}_{{ij}}$ $(1\leqslant i\lt i+2\leqslant j\leqslant N)$ and subsequently substitute Equation (22) into the resulting solution. Then, we obtain

$$\begin{eqnarray}{{\boldsymbol{P}}}_{{ij}} & = & {\tilde{{\boldsymbol{P}}}}_{{ij}}\equiv {m}_{i}{m}_{j}\\ & & \times (\displaystyle \sum _{k=i}^{j-1}\displaystyle \frac{1}{{m}_{k}{m}_{k+1}{| {{\boldsymbol{Q}}}_{k,k+1}| }^{2}}{{\boldsymbol{P}}}_{k,k+1}K{({{\boldsymbol{Q}}}_{k,k+1})}^{\top })K({{\boldsymbol{Q}}}_{{ij}})\\ & = & {m}_{i}{m}_{j}\\ & & \times [\begin{array}{c}{Q}_{{ij}[1]}{v}_{[1]}+{Q}_{{ij}[2]}{v}_{[2]}+{Q}_{{ij}[3]}{v}_{[3]}+{Q}_{{ij}[4]}{v}_{[4]}\\ {Q}_{{ij}[1]}{v}_{[2]}-{Q}_{{ij}[2]}{v}_{[1]}-{Q}_{{ij}[3]}{v}_{[4]}+{Q}_{{ij}[4]}{v}_{[3]}\\ {Q}_{{ij}[1]}{v}_{[3]}+{Q}_{{ij}[2]}{v}_{[4]}-{Q}_{{ij}[3]}{v}_{[1]}-{Q}_{{ij}[4]}{v}_{[2]}\\ -{Q}_{{ij}[1]}{v}_{[4]}+{Q}_{{ij}[2]}{v}_{[3]}-{Q}_{{ij}[3]}{v}_{[2]}+{Q}_{{ij}[4]}{v}_{[1]}\end{array}],\\ & & 1\leqslant i\lt i+2\leqslant j\leqslant N,\end{eqnarray} \tag{ 26 }$$

where

$$\begin{eqnarray*}{v}_{[1]} & \equiv & \displaystyle \sum _{k=i}^{j-1}({P}_{k,k+1[1]}{Q}_{k,k+1[1]}-{P}_{k,k+1[2]}{Q}_{k,k+1[2]}\\ & & -{P}_{k,k+1[3]}{Q}_{k,k+1[3]}+{P}_{k,k+1[4]}{Q}_{k,k+1[4]})/\\ & & ({m}_{k}{m}_{k+1}{| {{\boldsymbol{Q}}}_{k,k+1}| }^{2}),\\ {v}_{[2]} & \equiv & \displaystyle \sum _{k=i}^{j-1}({P}_{k,k+1[1]}{Q}_{k,k+1[2]}+{P}_{k,k+1[2]}{Q}_{k,k+1[1]}\\ & & -{P}_{k,k+1[3]}{Q}_{k,k+1[4]}-{P}_{k,k+1[4]}{Q}_{k,k+1[3]})/\\ & & ({m}_{k}{m}_{k+1}{| {{\boldsymbol{Q}}}_{k,k+1}| }^{2}),\\ {v}_{[3]} & \equiv & \displaystyle \sum _{k=i}^{j-1}({P}_{k,k+1[1]}{Q}_{k,k+1[3]}+{P}_{k,k+1[2]}{Q}_{k,k+1[4]}\\ & & +{P}_{k,k+1[3]}{Q}_{k,k+1[1]}+{P}_{k,k+1[4]}{Q}_{k,k+1[2]})/\\ & & ({m}_{k}{m}_{k+1}{| {{\boldsymbol{Q}}}_{k,k+1}| }^{2}),\\ {v}_{[4]} & \equiv & \displaystyle \sum _{k=i}^{j-1}({P}_{k,k+1[1]}{Q}_{k,k+1[4]}-{P}_{k,k+1[2]}{Q}_{k,k+1[3]}\\ & & +{P}_{k,k+1[3]}{Q}_{k,k+1[2]}-{P}_{k,k+1[4]}{Q}_{k,k+1[1]})/\\ & & ({m}_{k}{m}_{k+1}{| {{\boldsymbol{Q}}}_{k,k+1}| }^{2}).\end{eqnarray*}$$

Using Equations (22) and (26), we can rewrite Equation (17) as

$$\begin{eqnarray} & & \displaystyle \sum _{i=1}^{j-2}{\tilde{{\boldsymbol{G}}}}_{{ij}}-\displaystyle \sum _{i=j+2}^{N}{\tilde{{\boldsymbol{G}}}}_{{ji}}=-{{\boldsymbol{G}}}_{j-1,j}\\ & & +{{\boldsymbol{G}}}_{j,j+1},2\leqslant j\leqslant N-1;\\ & & \displaystyle \sum _{i=1}^{N-2}{\tilde{{\boldsymbol{G}}}}_{i,N}=-{{\boldsymbol{G}}}_{N-1,N},\end{eqnarray} \tag{ 27 }$$

where

$$\begin{eqnarray*}{\tilde{{\boldsymbol{G}}}}_{{ij}} & \equiv & \displaystyle \frac{1}{2{| {\tilde{{\boldsymbol{Q}}}}_{{ij}}| }^{2}}(\displaystyle \frac{d}{{dt}}{\tilde{{\boldsymbol{P}}}}_{{ij}}-(\displaystyle \frac{m}{4{m}_{i}{m}_{j}}{| {\tilde{{\boldsymbol{P}}}}_{{ij}}| }^{2}-2{m}_{i}{m}_{j})\displaystyle \frac{{\tilde{{\boldsymbol{Q}}}}_{{ij}}}{{| {\tilde{{\boldsymbol{Q}}}}_{{ij}}| }^{4}})\\ & & \quad {{\boldsymbol{K}}}^{\top }({\tilde{{\boldsymbol{Q}}}}_{{ij}}),1\leqslant i\lt i+2\leqslant j\leqslant N,\end{eqnarray*}$$

and ${{\boldsymbol{G}}}_{j-1,j}$, ${{\boldsymbol{G}}}_{j,j+1}$, and ${{\boldsymbol{G}}}_{N-1,N}$ are defined by Equation (14(b)). Note that Equation (27) is described by only $8(N-1)$ variables, ${Q}_{k,k+1[1]}$, ${Q}_{k,k+1[2]}$, ${Q}_{k,k+1[3]}$, ${Q}_{k,k+1[4]}$, ${P}_{k,k+1[1]}$, ${P}_{k,k+1[2]}$, ${P}_{k,k+1[3]}$, and ${P}_{k,k+1[4]}$, related to the chained vectors ${{\boldsymbol{q}}}_{k,k+1}$ and ${{\boldsymbol{p}}}_{k,k+1}$ $(1\leqslant k\leqslant N-1)$. In addition, the $2(N-1)$ vectors ${{\boldsymbol{Q}}}_{k,k+1}$ and ${{\boldsymbol{P}}}_{k,k+1}$ $(1\leqslant k\leqslant N-1)$ satisfy Equation (13(a)); namely,

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dt}}{{\boldsymbol{Q}}}_{k,k+1}=\displaystyle \frac{m}{4{m}_{k}{m}_{k+1}}\displaystyle \frac{{{\boldsymbol{P}}}_{k,k+1}}{{| {{\boldsymbol{Q}}}_{k,k+1}| }^{2}},\quad 1\leqslant k\leqslant N-1.\end{eqnarray} \tag{ 28 }$$

We call the system composed of Equations (27) and (28) the GNBP with chain variables. We clarify that the GNBP with chain variables describes the same motion as the RGNBP composed of Equations (13(a)), (15), and (17) in the following lemma.

Lemma 1. Suppose

• i.  
  The $2(N-1)$ vectors ${{\boldsymbol{Q}}}_{k,k+1}$ and ${{\boldsymbol{P}}}_{k,k+1}$ $(1\leqslant k\leqslant N-1)$ are the solutions of the GNBP with chain variables composed of Equations (27) and (28).
• ii.  
  The $(N-1)(N-2)$ vectors ${{\boldsymbol{Q}}}_{{ij}}$ and ${{\boldsymbol{P}}}_{{ij}}$ $(1\leqslant i\lt i+2\leqslant j\leqslant N)$ are given by Equations (22) and (26).

Then, the $N(N-1)$ vectors ${{\boldsymbol{Q}}}_{{ij}}$ and ${{\boldsymbol{P}}}_{{ij}}$ $(1\leqslant i\lt j\leqslant N)$ satisfy the RGNBP composed of Equations (13(a)), (15), and (17).

Proof. (a) Derivation of Equation (13(a)).

The vectors ${{\boldsymbol{Q}}}_{{ij}}$ and ${{\boldsymbol{P}}}_{{ij}}$ $(1\leqslant i\lt i+2\leqslant j\leqslant N)$ defined by Equations (22) and (26) satisfy Equation (25). In addition, using Equations (28), (25) is rewritten as

$$\begin{eqnarray}&&\displaystyle \frac{m}{4{m}_{i}{m}_{j}}\displaystyle \frac{{{\boldsymbol{P}}}_{{ij}}{\boldsymbol{K}}{({{\boldsymbol{Q}}}_{{ij}})}^{\top }}{{| {{\boldsymbol{Q}}}_{{ij}}| }^{2}}=\displaystyle \sum _{k=i}^{j-1}\displaystyle \frac{d}{{dt}}{{\boldsymbol{Q}}}_{k,k+1}{\boldsymbol{K}}{({{\boldsymbol{Q}}}_{k,k+1})}^{\top },\\ &&1\leqslant i\lt i+2\leqslant j\leqslant N.\end{eqnarray} \tag{ 29 }$$

Similarly, the vectors ${{\boldsymbol{Q}}}_{{ij}}$ $(1\leqslant i\lt i+2\leqslant j\leqslant N)$ in Equation (22) satisfy Equation (23), so the time derivative of Equation (23) is also fulfilled. The derivative is written in the following form:

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dt}}{{\boldsymbol{Q}}}_{{ij}}{\boldsymbol{K}}{({{\boldsymbol{Q}}}_{{ij}})}^{\top }=\displaystyle \sum _{k=i}^{j-1}\displaystyle \frac{d}{{dt}}{{\boldsymbol{Q}}}_{k,k+1}{\boldsymbol{K}}{({{\boldsymbol{Q}}}_{k,k+1})}^{\top },\\ &&1\leqslant i\lt i+2\leqslant j\leqslant N.\end{eqnarray} \tag{ 30 }$$

Because the right-hand side of Equation (29) coincides with that of Equation (30), the left-hand sides of Equations (29) and (30) are equal. Accordingly, Equation (13(a)) in the RGNBP is obtained.

(b) Derivation of Equation (15).

We have already shown that the vectors ${{\boldsymbol{Q}}}_{{ij}}$ $(1\leqslant i\lt i+2\leqslant j\leqslant N)$ in Equation (22) fulfill Equation (23). In addition, substituting Equation (23) into the right-hand side of ${{\boldsymbol{\Phi }}}_{1{jk}}({\boldsymbol{Q}})$ in Equation (15) yields $0$. Accordingly, the vectors ${{\boldsymbol{Q}}}_{{ij}}$ $(1\leqslant i\lt i+2\leqslant j\leqslant N)$ in Equation (22) satisfy Equation (15).

(c) Derivation of Equation (17).

Substituting Equations (22) and (26) into Equation (17) yields Equation (27). Namely, Equation (27) is equal to Equation (17) under the conditions in Equations (22) and (26). Thus, the $N(N-1)$ vectors ${{\boldsymbol{Q}}}_{{ij}}$ and ${{\boldsymbol{P}}}_{{ij}}$ $(1\leqslant i\lt j\leqslant N)$ satisfy Equation (17).□

For large N, the number of dependent variables ${P}_{i,i+1[1]}$, ${P}_{i,i+1[2]}$, ${P}_{i,i+1[3]}$, ${P}_{i,i+1[4]}$, ${Q}_{i,i+1[1]}$, ${Q}_{i,i+1[2]}$, ${Q}_{i,i+1[3]}$, and ${Q}_{i,i+1[4]}$ in the GNBP with chain variables, $8(N-1)$, is remarkably smaller than $4N(N-1)$, which is the number of dependent variables in the RGNBP composed of Equations (13(a)), (15), and (17). Therefore, the computational cost of the GNBP with chain variables is much lower than that of the RGNBP.

We apply the d'Alembert-type scheme (Betsch 2005) to the RGNBP composed of Equations ((13a)), (15), and (17). Then, the discrete-time system, which approximately describes the motion of the RGNBP in a typical time interval ${I}^{(n)}=[{t}^{(n)},{t}^{(n+1)}]$ with a corresponding time step ${\rm{\Delta }}{t}^{(n,n+1)}\equiv {t}^{(n+1)}-{t}^{(n)}$, is given as follows:

$$\begin{eqnarray}&&{{\boldsymbol{F}}}_{{ij}}^{(n+1)}=0,\quad 1\leqslant i\lt j\leqslant N,\end{eqnarray} \tag{ 31a }$$

$$\begin{eqnarray}&&\displaystyle \sum _{i=1}^{j-1}{{\boldsymbol{G}}}_{{ij}}^{(n+1)}-\displaystyle \sum _{i=j+1}^{N}{{\boldsymbol{G}}}_{{ji}}^{(n+1)}=0,2\leqslant j\leqslant N,\end{eqnarray} \tag{ 31b }$$

$$\begin{eqnarray}&&{{\boldsymbol{\Phi }}}_{1{ij}}({{\boldsymbol{Q}}}^{(n+1)})=0,\ 2\leqslant i\lt j\leqslant N,\end{eqnarray} \tag{ 31c }$$

where ${{\boldsymbol{Q}}}_{{ij}}^{(l)}=({Q}_{{ij}[1]}^{(l)},{Q}_{{ij}[2]}^{(l)},{Q}_{{ij}[3]}^{(l)},{Q}_{{ij}[4]}^{(l)})$ and ${{\boldsymbol{P}}}_{{ij}}^{(l)}=({P}_{{ij}[1]}^{(l)},{P}_{{ij}[2]}^{(l)},{P}_{{ij}[3]}^{(l)},{P}_{{ij}[4]}^{(l)})$ at time ${t}^{(l)}$ $(l=n,n+1)$ for $1\leqslant i\lt j\leqslant N$. We define the midpoint value ${(\bullet )}^{(n+1/2)}\;\equiv [{(\bullet )}^{(n+1)}+{(\bullet )}^{(n)}]/2$ of the function $(\bullet )(t)$, and

$$\begin{eqnarray}{{\boldsymbol{F}}}_{{ij}}^{(n+1)} & \equiv & \displaystyle \frac{{{\boldsymbol{Q}}}_{{ij}}^{(n+1)}-{{\boldsymbol{Q}}}_{{ij}}^{(n)}}{{\rm{\Delta }}{t}^{(n,n+1)}}\\ & & -\displaystyle \frac{m}{8{m}_{i}{m}_{j}}\displaystyle \frac{{| {{\boldsymbol{Q}}}_{{ij}}^{(n+1)}| }^{2}+{| {{\boldsymbol{Q}}}_{{ij}}^{(n)}| }^{2}}{{| {{\boldsymbol{Q}}}_{{ij}}^{(n+1)}| }^{2}{| {{\boldsymbol{Q}}}_{{ij}}^{(n)}| }^{2}}{{\boldsymbol{P}}}_{{ij}}^{(n+1/2)},\\ & & 1\leqslant i\lt j\leqslant N,\end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray}{{\boldsymbol{G}}}_{{ij}}^{(n+1)} & \equiv & \displaystyle \frac{1}{2{| {{\boldsymbol{Q}}}_{{ij}}^{(n+1/2)}| }^{2}}(\displaystyle \frac{{{\boldsymbol{P}}}_{{ij}}^{(n+1)}-{{\boldsymbol{P}}}_{{ij}}^{(n)}}{{\rm{\Delta }}{t}^{(n,n+1)}}\\ & & -\displaystyle \frac{1}{{| {{\boldsymbol{Q}}}_{{ij}}^{(n+1)}| }^{2}{| {{\boldsymbol{Q}}}_{{ij}}^{(n)}| }^{2}}\\ & & (\displaystyle \frac{m}{8{m}_{i}{m}_{j}}({| {{\boldsymbol{P}}}_{{ij}}^{(n+1)}| }^{2}+{| {{\boldsymbol{P}}}_{{ij}}^{(n)}| }^{2})-2{m}_{i}{m}_{j}){{\boldsymbol{Q}}}_{{ij}}^{(n+1/2)})\\ & & {\boldsymbol{K}}({{\boldsymbol{Q}}}_{{ij}}^{(n+1/2)}){}^{\top },\;1\leqslant i\lt j\leqslant N.\end{eqnarray} \tag{ 33 }$$

${{\boldsymbol{Q}}}^{(n)}\in {\mathcal{Q}}$ and ${{\boldsymbol{P}}}^{(n)}\in {{\mathbb{R}}}^{2N(N-1)}$ are assumed to be quantities at time node ${t}^{(n)}$, where ${\mathcal{Q}}$ was already defined by Equation (19). In the following, we call this system described by Equation (31) the discrete-time redundant general N-body problem (d-RGNBP). It can be used to calculate the unknown ${{\boldsymbol{Q}}}^{(n+1)}$ and $\begin{array}{ccccc}{{\boldsymbol{P}}}^{(n+1)}=({{\boldsymbol{P}}}_{1,2}^{(n+1)}, & \cdots & {{\boldsymbol{P}}}_{1,N}^{(n+1)},| {{\boldsymbol{P}}}_{2,3}^{(n+1)}, & \cdots & {{\boldsymbol{P}}}_{2,N}^{(n+1)},\end{array}$ $| \cdots | {{\boldsymbol{P}}}_{N-1,N}^{(n+1)})\in {{\mathbb{R}}}^{2N(N-1)}$. Note that ${{\boldsymbol{Q}}}^{(n+1)}\in {\mathcal{Q}}$ because of Equation ((31c)).

The d-RGNBP precisely preserves all the conserved quantities except the angular momentum, as demonstrated in the following theorem.

Theorem 1. (Conserved quantities of d-RGNBP) The d-RGNBP in Equation (31) preserves the following three conserved quantities exactly:

• 1.  
  The value of the Hamiltonian ${h}_{\mathrm{rel}}$ defined by Equation (9),
• 2.  
  The linear momentum ${\boldsymbol{l}}\equiv {\displaystyle \sum }_{i=1}^{N}{{\boldsymbol{p}}}_{i}^{\prime }=0$,
• 3.  
  The position of the center of mass ${\boldsymbol{c}}\equiv {\displaystyle \sum }_{i=1}^{N}{m}_{i}{{\boldsymbol{q}}}_{i}^{\prime }=0$.

Proof. 

• 1.  
  Using the multipliers ${\bar{{\boldsymbol{\Lambda }}}}_{{ij}}=({{\rm{\Lambda }}}_{{ij}[1]},{{\rm{\Lambda }}}_{{ij}[2]},{{\rm{\Lambda }}}_{{ij}[3]},0)$ $(2\leqslant i\lt j\leqslant N)$, the d-RGNBP in Equation (31) is rewritten as

  $$\begin{eqnarray}&&{{\boldsymbol{F}}}_{{ij}}^{(n+1)}=0,\quad 1\leqslant i\lt j\leqslant N,\end{eqnarray} \tag{ 34a }$$

  $$\begin{eqnarray}&&{{\boldsymbol{G}}}_{1j}^{(n+1)}-\displaystyle \sum _{k=2}^{j-1}{\bar{{\boldsymbol{\Lambda }}}}_{{kj}}+\displaystyle \sum _{k=j+1}^{N}{\bar{{\boldsymbol{\Lambda }}}}_{{jk}}=0,\quad j=2,\cdots ,N,\end{eqnarray} \tag{ 34b }$$

  $$\begin{eqnarray}&&{{\boldsymbol{G}}}_{{ij}}^{(n+1)}+{\bar{{\boldsymbol{\Lambda }}}}_{{ij}}=0,\quad 2\leqslant i\lt j\leqslant N,\end{eqnarray} \tag{ 34c }$$

  $$\begin{eqnarray}&&{{\boldsymbol{\Phi }}}_{1{ij}}({{\boldsymbol{Q}}}_{{ij}}^{(n+1)})=0,\quad 2\leqslant i\lt j\leqslant N.\end{eqnarray} \tag{ 34d }$$

  Using Equations (34(c)), (34(b)) is expressed as Equation (31(b)). We set the following functions:

  $$\begin{eqnarray*}&&{H}_{1j}({{\boldsymbol{P}}}_{1j},{{\boldsymbol{Q}}}_{1j})\equiv \displaystyle \frac{m}{8{m}_{1}{m}_{j}}\displaystyle \frac{{| {{\boldsymbol{P}}}_{1j}| }^{2}}{{| {{\boldsymbol{Q}}}_{1j}| }^{2}}-\displaystyle \frac{{m}_{1}{m}_{j}}{{| {{\boldsymbol{Q}}}_{1j}| }^{2}}\\ &&+{{\boldsymbol{Q}}}_{1j}{\boldsymbol{K}}{({{\boldsymbol{Q}}}_{1j})}^{\top }\cdot (\displaystyle \sum _{i=2}^{j-1}{\bar{{\boldsymbol{\Lambda }}}}_{{ij}}-\displaystyle \sum _{i=j+1}^{N}{\bar{{\boldsymbol{\Lambda }}}}_{{ji}}),2\leqslant j\leqslant N,\\ &&{H}_{{ij}}({{\boldsymbol{P}}}_{{ij}},{{\boldsymbol{Q}}}_{{ij}})\equiv \displaystyle \frac{m}{8{m}_{i}{m}_{j}}\displaystyle \frac{{| {{\boldsymbol{P}}}_{{ij}}| }^{2}}{{| {{\boldsymbol{Q}}}_{{ij}}| }^{2}}-\displaystyle \frac{{m}_{i}{m}_{j}}{{| {{\boldsymbol{Q}}}_{{ij}}| }^{2}}\\ &&+{{\boldsymbol{Q}}}_{{ij}}{\boldsymbol{K}}{({{\boldsymbol{Q}}}_{{ij}})}^{\top }\cdot {\bar{{\boldsymbol{\Lambda }}}}_{{ij}},\quad 2\leqslant i\lt j\leqslant N.\end{eqnarray*}$$

  Because of Equation (15) and

  $$\begin{eqnarray*}&&\displaystyle \sum _{i=1}^{N-1}\displaystyle \sum _{j=i+1}^{N}{H}_{{ij}}({{\boldsymbol{P}}}_{{ij}},{{\boldsymbol{Q}}}_{{ij}})\\ &&\quad =\displaystyle \sum _{i=1}^{N-1}\displaystyle \sum _{j=i+1}^{N}(\displaystyle \frac{m}{8{m}_{i}{m}_{j}}\displaystyle \frac{{| {{\boldsymbol{P}}}_{{ij}}| }^{2}}{{| {{\boldsymbol{Q}}}_{{ij}}| }^{2}}-\displaystyle \frac{{m}_{i}{m}_{j}}{{| {{\boldsymbol{Q}}}_{{ij}}| }^{2}})\\ &&\qquad +\displaystyle \sum _{j=2}^{N-1}\displaystyle \sum _{k=j+1}^{N}{{\boldsymbol{\Phi }}}_{1{jk}}({\boldsymbol{Q}})\cdot {\bar{{\boldsymbol{\Lambda }}}}_{{jk}},\end{eqnarray*}$$

  the value of the right-hand side coincides with ${h}_{{\rm{LC}}}$ in Equation (18), which equals the value of ${H}_{{\rm{LC}}}$ defined by Equation (16). For $2\leqslant j\leqslant N$, the scalar multiplication of Equation (34(a)) for i = 1 and Equation (34(b)) by ${{\boldsymbol{P}}}_{1j}^{(n+1)}-{{\boldsymbol{P}}}_{1j}^{(n)}$ and $-({{\boldsymbol{Q}}}_{1j}^{(n+1)}-{{\boldsymbol{Q}}}_{1j}^{(n)})$, respectively, and the subsequent addition of the two equations yield

  $$\begin{eqnarray} & & \displaystyle \frac{m}{8{m}_{1}{m}_{j}}\displaystyle \frac{{| {{\boldsymbol{P}}}_{1j}^{(n+1)}| }^{2}}{{| {{\boldsymbol{Q}}}_{1j}^{(n+1)}| }^{2}}\\ & & -\displaystyle \frac{m}{8{m}_{1}{m}_{j}}\displaystyle \frac{{| {{\boldsymbol{P}}}_{1j}^{(n)}| }^{2}}{{| {{\boldsymbol{Q}}}_{1j}^{(n)}| }^{2}}-(\displaystyle \frac{{m}_{1}{m}_{j}}{{| {{\boldsymbol{Q}}}_{1j}^{(n+1)}| }^{2}}-\displaystyle \frac{{m}_{1}{m}_{j}}{{| {{\boldsymbol{Q}}}_{1j}^{(n)}| }^{2}})\\ & & +({{\boldsymbol{Q}}}_{1j}^{(n+1)}{\boldsymbol{K}}{({{\boldsymbol{Q}}}_{1j}^{(n+1)})}^{\top }-{{\boldsymbol{Q}}}_{1j}^{(n)}{\boldsymbol{K}}{({{\boldsymbol{Q}}}_{1j}^{(n)})}^{\top })\\ & \cdot & (\displaystyle \sum _{i=2}^{j-1}{\bar{{\boldsymbol{\Lambda }}}}_{{ij}}-\displaystyle \sum _{i=j+1}^{N}{\bar{{\boldsymbol{\Lambda }}}}_{{ji}})=0,2\leqslant j\leqslant N.\end{eqnarray} \tag{ 35 }$$

  Equation (35) represents the conservation of ${H}_{1j}({{\boldsymbol{P}}}_{1j},{{\boldsymbol{Q}}}_{1j})$ because it equals ${H}_{1j}({{\boldsymbol{P}}}_{1j}^{(n+1)},{{\boldsymbol{Q}}}_{1j}^{(n+1)})\;-{H}_{1j}({{\boldsymbol{P}}}_{1j}^{(n)},{{\boldsymbol{Q}}}_{1j}^{(n)})=0$. In addition, for $2\leqslant i\lt j\leqslant N$, scalar multiplication of Equation (34(a)) for $i\geqslant 2$ and Equation (34(c)) by ${{\boldsymbol{P}}}_{{ij}}^{(n+1)}-{{\boldsymbol{P}}}_{{ij}}^{(n)}$ and $-({{\boldsymbol{Q}}}_{{ij}}^{(n+1)}-{{\boldsymbol{Q}}}_{{ij}}^{(n)})$, respectively, and subsequent addition of the two equations yield

  $$\begin{eqnarray} & & \displaystyle \frac{m}{8{m}_{i}{m}_{j}}\displaystyle \frac{{| {{\boldsymbol{P}}}_{{ij}}^{(n+1)}| }^{2}}{{| {{\boldsymbol{Q}}}_{{ij}}^{(n+1)}| }^{2}}-\displaystyle \frac{m}{8{m}_{i}{m}_{j}}\displaystyle \frac{{| {{\boldsymbol{P}}}_{{ij}}^{(n)}| }^{2}}{{| {{\boldsymbol{Q}}}_{{ij}}^{(n)}| }^{2}}\\ & & -(\displaystyle \frac{{m}_{i}{m}_{j}}{{| {{\boldsymbol{Q}}}_{{ij}}^{(n+1)}| }^{2}}-\displaystyle \frac{{m}_{i}{m}_{j}}{{| {{\boldsymbol{Q}}}_{{ij}}^{(n)}| }^{2}})\\ & & +({{\boldsymbol{Q}}}_{{ij}}^{(n+1)}{\boldsymbol{K}}{({{\boldsymbol{Q}}}_{{ij}}^{(n+1)})}^{\top }-{{\boldsymbol{Q}}}_{{ij}}^{(n)}{\boldsymbol{K}}{({{\boldsymbol{Q}}}_{{ij}}^{(n)})}^{\top })\\ & & \cdot {\bar{{\boldsymbol{\Lambda }}}}_{{ij}}=0,2\leqslant i\lt j\leqslant N.\end{eqnarray} \tag{ 36 }$$

  Equation (36) describes the conservation of ${H}_{{ij}}({{\boldsymbol{P}}}_{{ij}},{{\boldsymbol{Q}}}_{{ij}})$ because it is equivalent to ${H}_{{ij}}({{\boldsymbol{P}}}_{{ij}}^{(n+1)},{{\boldsymbol{Q}}}_{{ij}}^{(n+1)})\;-{H}_{{ij}}({{\boldsymbol{P}}}_{{ij}}^{(n)},{{\boldsymbol{Q}}}_{{ij}}^{(n)})=0$. Accordingly, the d-RGNBP in Equation (31) preserves the value of ${h}_{{\rm{LC}}}$ defined by Equation (18). Because the value of ${h}_{{\rm{LC}}}$ is that of ${H}_{{\rm{LC}}}$ in Equation (16), and the value of ${h}_{{\rm{LC}}}$ equals ${h}_{\mathrm{rel}}$ in Equation (9), as described in Section 2.2.2, the d-RGNBP in Equation (31) conserves the value of ${H}_{\mathrm{rel}}$ defined by Equation (9).
• 2.  
  The linear momentum ${{\boldsymbol{l}}}^{(n+1)}$ in the barycentric frame at time node ${t}^{(n+1)}$ is given by

  $$\begin{eqnarray}&&{{\boldsymbol{l}}}^{(n+1)}=\displaystyle \sum _{i=1}^{N}\;{{\boldsymbol{p}}}_{i}^{\prime (n+1)}.\end{eqnarray} \tag{ 37 }$$

  Substituting Equation (6) into Equation (37), we see

  $$\begin{eqnarray*}&&{{\boldsymbol{l}}}^{(n+1)}=0.\end{eqnarray*}$$

  Therefore, the d-RGNBP in Equation (31) preserves the linear momentum ${\boldsymbol{l}}$, which is identically zero.
• 3.  
  The position of the center of mass ${{\boldsymbol{c}}}^{(n+1)}$ in the barycentric frame at time node ${t}^{(n+1)}$ is expressed by

  $$\begin{eqnarray}&&{{\boldsymbol{c}}}^{(n+1)}=\displaystyle \sum _{i=1}^{N}{m}_{i}{{\boldsymbol{q}}}_{i}^{\prime (n+1)}.\end{eqnarray} \tag{ 38 }$$

  We substitute Equation (6) into (38). Then, we can check that

  $$\begin{eqnarray*}&&{{\boldsymbol{c}}}^{(n+1)}=0.\end{eqnarray*}$$

As a result, the d-RGNBP in Equation (31) maintains the vector value of the position of the center of motion ${\boldsymbol{c}}$ at zero.□

For large N, the number of variables ${P}_{{ij}[1]}^{(n+1)}$, ${P}_{{ij}[2]}^{(n+1)}$, ${P}_{{ij}[3]}^{(n+1)}$, ${P}_{{ij}[4]}^{(n+1)}$, ${Q}_{{ij}[1]}^{(n+1)}$, ${Q}_{{ij}[2]}^{(n+1)}$, ${Q}_{{ij}[3]}^{(n+1)}$, and ${Q}_{{ij}[4]}^{(n+1)}$ in the d-RGNBP in Equation (31), $4N(N-1)$, is remarkably larger than $6N$, which is the number of dependent variables $p{^{\prime} }_{i[1]}$, $p{^{\prime} }_{i[2]}$, $p{^{\prime} }_{i[3]}$, ${q}_{i[1]}^{\prime }$, ${q}_{i[2]}^{\prime }$, and ${q}_{i[3]}^{\prime }$ in the GNBP in Equation (2). Thus, the computational cost of the d-RGNBP in Equation (31) is very high.

In this section, to reduce the high computational cost of the d-RGNBP, we rewrite it using only ${{\boldsymbol{Q}}}_{k,k+1}^{(n)}$, ${{\boldsymbol{Q}}}_{k,k+1}^{(n+1)}$, ${{\boldsymbol{P}}}_{k,k+1}^{(n)}$, and ${{\boldsymbol{P}}}_{k,k+1}^{(n+1)}$ $(1\leqslant k\leqslant N-1)$ related to chained vectors. First, we show that the vectors ${{\boldsymbol{Q}}}_{{ij}}^{(l)}$ $(1\leqslant i\lt i+2\leqslant j\leqslant N,l=n,n+1)$ can be expressed as functions of ${{\boldsymbol{Q}}}_{k,k+1}^{(l)}$ $(1\leqslant k\leqslant N-1,l=n,n+1)$ related to the chained position vectors ${{\boldsymbol{q}}}_{k,k+1}^{(l)}$ $(1\leqslant k\leqslant N-1,l=n,n+1)$.

Because ${{\boldsymbol{Q}}}^{(l)}\in {\mathcal{Q}}$ $(l=n,n+1)$, where ${\mathcal{Q}}$ is a $4(N-1)$-dimensional manifold defined by Equation (19), we can assume that the $(N-1)$ vectors ${{\boldsymbol{Q}}}_{k,k+1}^{(l)}$ $(1\leqslant k\leqslant N-1,l=n,n+1)$, which correspond to $4(N-1)$ variables, constitute a basis for the manifold ${\mathcal{Q}}$. According to Equation (10), from these base vectors, the chained position vectors are given by

$$\begin{eqnarray}&&{{\boldsymbol{q}}}_{k,k+1}^{(l)}\equiv {[\begin{array}{c}{q}_{k,k+1[1]}^{(l)}\\ {q}_{k,k+1[2]}^{(l)}\\ {q}_{k,k+1[3]}^{(l)}\end{array}]}^{\top }\\ &&=\;{[\begin{array}{c}{({Q}_{k,k+1[1]}^{(l)})}^{2}-\;{({Q}_{k,k+1[2]}^{(l)})}^{2}\\ \quad -{({Q}_{k,k+1[3]}^{(l)})}^{2}+{({Q}_{k,k+1[4]}^{(l)})}^{2}\\ 2({Q}_{k,k+1[1]}^{(l)}{Q}_{k,k+1[2]}^{(l)}-{Q}_{k,k+1[3]}^{(l)}{Q}_{k,k+1[4]}^{(l)})\\ 2({Q}_{k,k+1[1]}^{(l)}{Q}_{k,k+1[3]}^{(l)}+{Q}_{k,k+1[2]}^{(l)}{Q}_{k,k+1[4]}^{(l)})\end{array}]}^{\top },\\ &&1\leqslant k\leqslant N-1,l=n,n+1.\end{eqnarray} \tag{ 39 }$$

In this basis, each vector ${{\boldsymbol{Q}}}_{{ij}}^{(l)}\in {\mathcal{Q}}$ $(1\leqslant i\lt i+2\leqslant j\leqslant N,l=n,n+1)$ related to a non-chained position vector ${{\boldsymbol{q}}}_{{ij}}^{(l)}$ satisfies

$$\begin{eqnarray}{{\boldsymbol{Q}}}_{{ij}}^{(l)}{\boldsymbol{K}}{({{\boldsymbol{Q}}}_{{ij}}^{(l)})}^{\top } & = & \displaystyle \sum _{k=i}^{j-1}{{\boldsymbol{Q}}}_{k,k+1}^{(l)}{\boldsymbol{K}}{({{\boldsymbol{Q}}}_{k,k+1}^{(l)})}^{\top },\\ & & \quad 1\leqslant i\lt i+2\leqslant j\leqslant N,l=n,n+1.\end{eqnarray} \tag{ 40 }$$

In addition, to avoid cancellation of significant digits, the solutions of Equation (40) for ${{\boldsymbol{Q}}}_{{ij}}^{(l)}$ are written as follows:

$$\begin{eqnarray}{{\boldsymbol{Q}}}_{{ij}}^{(l)} & = & {\tilde{{\boldsymbol{Q}}}}_{{ij}}^{(l)}\equiv ({\tilde{Q}}_{{ij}[1]}^{(l)},{\tilde{Q}}_{{ij}[2]}^{(l)},{\tilde{Q}}_{{ij}[3]}^{(l)},{\tilde{Q}}_{{ij}[4]}^{(l)})\\ & \equiv & \{\begin{array}{cc}({\tilde{u}}_{{ij}[1]}^{(l)},{\tilde{u}}_{{ij}[2]}^{(l)},{\tilde{u}}_{{ij}[3]}^{(l)},{\tilde{u}}_{{ij}[4]}^{(l)}) & \begin{array}{l}\mathrm{if}\ \displaystyle \sum _{k=i}^{j-1}{q}_{k,k+1[1]}^{(l)}\geqslant 0,\\ l=n,n+1,\end{array}\\ ({\tilde{u}}_{{ij}[2]}^{(l)},{\tilde{u}}_{{ij}[1]}^{(l)},{\tilde{u}}_{{ij}[4]}^{(l)},{\tilde{u}}_{{ij}[3]}^{(l)}) & \begin{array}{l}\mathrm{if}\ \displaystyle \sum _{k=i}^{j-1}{q}_{k,k+1[1]}^{(l)}\lt 0,\\ l=n,n+1,\end{array}\end{array}\end{eqnarray} \tag{ 41 }$$

where

$$\begin{eqnarray}{\tilde{u}}_{{ij}[1]}^{(l)} & = & {(\displaystyle \frac{| \displaystyle \sum _{k=i}^{j-1}{{\boldsymbol{q}}}_{k,k+1}^{(l)}| +| \displaystyle \sum _{k=i}^{j-1}{q}_{k,k+1[1]}^{(l)}| }{2})}^{1/2},\\ & & \quad {\tilde{u}}_{{ij}[2]}^{(l)}=\displaystyle \frac{\displaystyle \sum _{k=i}^{j-1}{q}_{k,k+1[2]}^{(l)}}{2{\tilde{u}}_{{ij}[1]}^{(l)}},\\ {\tilde{u}}_{{ij}[3]}^{(l)} & = & \displaystyle \frac{\displaystyle \sum _{k=i}^{j-1}{q}_{k,k+1[3]}^{(l)}}{2{\tilde{u}}_{{ij}[1]}^{(l)}},\quad {\tilde{u}}_{{ij}[4]}^{(l)}=0,\\ & & \quad 1\leqslant i\lt i+2\leqslant j\leqslant N,l=n,n+1,\end{eqnarray} \tag{ 42 }$$

and ${{\boldsymbol{q}}}_{k,k+1}^{(l)}=({q}_{k,k+1[1]}^{(l)},{q}_{k,k+1[2]}^{(l)},{q}_{k,k+1[3]}^{(l)})$ in Equation (42) has been already defined by Equation (39).

Next, we write the vectors ${{\boldsymbol{P}}}_{{ij}}^{(n+1/2)}$ $(1\leqslant i\lt i+2\leqslant j\leqslant N)$ related to the non-chained momentum vectors ${{\boldsymbol{p}}}_{{ij}}^{(n+1/2)}$ as functions of ${{\boldsymbol{P}}}_{k,k+1}^{(n+1/2)}$, ${{\boldsymbol{Q}}}_{k,k+1}^{(n)}$, and ${{\boldsymbol{Q}}}_{k,k+1}^{(n+1)}$ $(1\leqslant k\leqslant N-1)$. Subtraction of Equation (40) with $l=n+1$ from that with l = n and subsequent division by $2{\rm{\Delta }}{t}^{(n,n+1)}=2({t}^{(n+1)}-{t}^{(n)})$ yield

$$\begin{eqnarray*} & & \displaystyle \frac{{{\boldsymbol{Q}}}_{{ij}}^{(n+1)}{\boldsymbol{K}}{({{\boldsymbol{Q}}}_{{ij}}^{(n+1)})}^{\top }-{{\boldsymbol{Q}}}_{{ij}}^{(n)}{\boldsymbol{K}}{({{\boldsymbol{Q}}}_{{ij}}^{(n)})}^{\top }}{2{\rm{\Delta }}{t}^{(n,n+1)}}\\ & = & \displaystyle \sum _{k=i}^{j-1}\displaystyle \frac{{{\boldsymbol{Q}}}_{k,k+1}^{(n+1)}{\boldsymbol{K}}{({{\boldsymbol{Q}}}_{k,k+1}^{(n+1)})}^{\top }-{{\boldsymbol{Q}}}_{k,k+1}^{(n)}{\boldsymbol{K}}{({{\boldsymbol{Q}}}_{k,k+1}^{(n)})}^{\top }}{2{\rm{\Delta }}{t}^{(n,n+1)}},\\ & & 1\leqslant i\lt i+2\leqslant j\leqslant N.\end{eqnarray*}$$

This leads to

$$\begin{eqnarray} & & \displaystyle \frac{{{\boldsymbol{Q}}}_{{ij}}^{(n+1)}-{{\boldsymbol{Q}}}_{{ij}}^{(n)}}{{\rm{\Delta }}{t}^{(n,n+1)}}{\boldsymbol{K}}{({{\boldsymbol{Q}}}_{{ij}}^{(n+1/2)})}^{\top }\\ & = & \displaystyle \sum _{k=i}^{j-1}\displaystyle \frac{{{\boldsymbol{Q}}}_{k,k+1}^{(n+1)}-{{\boldsymbol{Q}}}_{k,k+1}^{(n)}}{{\rm{\Delta }}{t}^{(n,n+1)}}{\boldsymbol{K}}{({{\boldsymbol{Q}}}_{k,k+1}^{(n+1/2)})}^{\top },\\ & & 1\leqslant i\lt i+2\leqslant j\leqslant N.\end{eqnarray} \tag{ 43 }$$

Using Equations (31(a)) and (41), Equation (43) is rewritten as

$$\begin{eqnarray}&&\displaystyle \frac{m}{8{m}_{i}{m}_{j}}\displaystyle \frac{{| {\tilde{{\boldsymbol{Q}}}}_{{ij}}^{(n+1)}| }^{2}+{| {\tilde{{\boldsymbol{Q}}}}_{{ij}}^{(n)}| }^{2}}{{| {\tilde{{\boldsymbol{Q}}}}_{{ij}}^{(n+1)}| }^{2}{| {\tilde{{\boldsymbol{Q}}}}_{{ij}}^{(n)}| }^{2}}{{\boldsymbol{P}}}_{{ij}}^{(n+1/2)}{\boldsymbol{K}}{({\tilde{{\boldsymbol{Q}}}}_{{ij}}^{(n+1/2)})}^{\top }\\ &&=\displaystyle \sum _{k=i}^{j-1}\displaystyle \frac{m}{8{m}_{k}{m}_{k+1}}\displaystyle \frac{{| {{\boldsymbol{Q}}}_{k,k+1}^{(n+1)}| }^{2}+{| {{\boldsymbol{Q}}}_{k,k+1}^{(n)}| }^{2}}{{| {{\boldsymbol{Q}}}_{k,k+1}^{(n+1)}| }^{2}{| {{\boldsymbol{Q}}}_{k,k+1}^{(n)}| }^{2}}{{\boldsymbol{P}}}_{k,k+1}^{(n+1/2)}{\boldsymbol{K}}{({{\boldsymbol{Q}}}_{k,k+1}^{(n+1/2)})}^{\top },\\ &&1\leqslant i\lt i+2\leqslant j\leqslant N.\end{eqnarray} \tag{ 44 }$$

We solve Equation (44) for ${{\boldsymbol{P}}}_{{ij}}^{(n+1/2)}$ $(1\leqslant i\lt i+2\leqslant j\leqslant N)$ and subsequently substitute Equation (41) into the resulting solution. Then, we have

$$\begin{eqnarray}&&{{\boldsymbol{P}}}_{{ij}}^{(n+1/2)}={\tilde{{\boldsymbol{P}}}}_{{ij}}^{(n+1/2)}\\ &&\equiv {m}_{i}{m}_{j}\displaystyle \frac{{| {\tilde{{\boldsymbol{Q}}}}_{{ij}}^{(n+1)}| }^{2}{| {\tilde{{\boldsymbol{Q}}}}_{{ij}}^{(n)}| }^{2}}{{| {\tilde{{\boldsymbol{Q}}}}_{{ij}}^{(n+1/2)}| }^{2}({| {\tilde{{\boldsymbol{Q}}}}_{{ij}}^{(n+1)}| }^{2}+{| {\tilde{{\boldsymbol{Q}}}}_{{ij}}^{(n)}| }^{2})}\\ &&(\displaystyle \sum _{k=i}^{j-1}\displaystyle \frac{1}{{m}_{k}{m}_{k+1}}\displaystyle \frac{{| {{\boldsymbol{Q}}}_{k,k+1}^{(n+1)}| }^{2}+{| {{\boldsymbol{Q}}}_{k,k+1}^{(n)}| }^{2}}{{| {{\boldsymbol{Q}}}_{k,k+1}^{(n+1)}| }^{2}{| {{\boldsymbol{Q}}}_{k,k+1}^{(n)}| }^{2}}{{\boldsymbol{P}}}_{k,k+1}^{(n+1/2)}{\boldsymbol{K}}{({{\boldsymbol{Q}}}_{k,k+1}^{(n+1/2)})}^{\top })\\ &&{\boldsymbol{K}}({\tilde{{\boldsymbol{Q}}}}_{{ij}}^{(n+1/2)}),1\leqslant i\lt i+2\leqslant j\leqslant N.\end{eqnarray} \tag{ 45 }$$

Because of the definition of the midpoint value ${(\bullet )}^{(n+1/2)}$ and Equation (45),

$$\begin{eqnarray}&&{{\boldsymbol{P}}}_{{ij}}^{(n+1)}={\tilde{{\boldsymbol{P}}}}_{{ij}}^{(n+1)}\equiv 2{\tilde{{\boldsymbol{P}}}}_{{ij}}^{(n+1/2)}-{\tilde{{\boldsymbol{P}}}}_{{ij}}^{(n)}\ (1\leqslant i\lt j\leqslant N).\end{eqnarray} \tag{ 46 }$$

Using Equations (41) and (45), we can rewrite Equation (31(b)) in the d-RGNBP as

$$\begin{eqnarray}&&\displaystyle \sum _{i=1}^{j-2}{\tilde{{\boldsymbol{G}}}}_{{ij}}^{(n+1)}-\displaystyle \sum _{i=j+2}^{N}{\tilde{{\boldsymbol{G}}}}_{{ji}}^{(n+1)}\\ &&=-{{\boldsymbol{G}}}_{j-1,j}^{(n+1)}+{{\boldsymbol{G}}}_{j,j+1}^{(n+1)},\ 2\leqslant j\leqslant N-1;\\ &&\quad \displaystyle \sum _{i=1}^{N-2}{\tilde{{\boldsymbol{G}}}}_{i,N}^{(n+1)}=-{{\boldsymbol{G}}}_{N-1,N}^{(n+1)},\end{eqnarray} \tag{ 47 }$$

where

$$\begin{eqnarray*}{\tilde{{\boldsymbol{G}}}}_{{ij}}^{(n+1)} & \equiv & \displaystyle \frac{1}{2{| {\tilde{{\boldsymbol{Q}}}}_{{ij}}^{(n+1/2)}| }^{2}}(\displaystyle \frac{{\tilde{{\boldsymbol{P}}}}_{{ij}}^{(n+1)}-{\tilde{{\boldsymbol{P}}}}_{{ij}}^{(n)}}{{\rm{\Delta }}{t}^{(n,n+1)}}-\displaystyle \frac{1}{{| {\tilde{{\boldsymbol{Q}}}}_{{ij}}^{(n+1)}| }^{2}{| {\tilde{{\boldsymbol{Q}}}}_{{ij}}^{(n)}| }^{2}}\\ & & \times (\displaystyle \frac{m}{8{m}_{i}{m}_{j}}({| {\tilde{{\boldsymbol{P}}}}_{{ij}}^{(n+1)}| }^{2}+{| {\tilde{{\boldsymbol{P}}}}_{{ij}}^{(n)}| }^{2})\\ & & -2{m}_{i}{m}_{j}){\tilde{{\boldsymbol{Q}}}}_{{ij}}^{(n+1/2)}){\boldsymbol{K}}{({\tilde{{\boldsymbol{Q}}}}_{{ij}}^{(n+1/2)})}^{\top },\\ & & 1\leqslant i\lt i+2\leqslant j\leqslant N,\end{eqnarray*}$$

and ${{\boldsymbol{G}}}_{j-1,j}^{(n+1)}$, ${{\boldsymbol{G}}}_{j,j+1}^{(n+1)}$, and ${{\boldsymbol{G}}}_{N-1,N}^{(n+1)}$ have the same form as Equation (33). Note that subsystem (47) is described by only $4(N-1)$ vectors, ${{\boldsymbol{Q}}}_{k,k+1}^{(n)}$, ${{\boldsymbol{Q}}}_{k,k+1}^{(n+1)}$, ${{\boldsymbol{P}}}_{k,k+1}^{(n)}$, and ${{\boldsymbol{P}}}_{k,k+1}^{(n+1)}$. In addition, these vectors satisfy Equation (31(a)); namely,

$$\begin{eqnarray}&&\displaystyle \frac{{{\boldsymbol{Q}}}_{k,k+1}^{(n+1)}-{{\boldsymbol{Q}}}_{k,k+1}^{(n)}}{{\rm{\Delta }}{t}^{(n,n+1)}}\\ &&\quad =\displaystyle \frac{m}{8{m}_{k}{m}_{k+1}}\displaystyle \frac{{| {{\boldsymbol{Q}}}_{k,k+1}^{(n+1)}| }^{2}+{| {{\boldsymbol{Q}}}_{k,k+1}^{(n)}| }^{2}}{{| {{\boldsymbol{Q}}}_{k,k+1}^{(n+1)}| }^{2}{| {{\boldsymbol{Q}}}_{k,k+1}^{(n)}| }^{2}}{{\boldsymbol{P}}}_{k,k+1}^{(n+1/2)},\\ &&1\leqslant k\leqslant N-1.\end{eqnarray} \tag{ 48 }$$

We call the discrete-time system composed of Equations (47) and (48) the discrete-time GNBP with chain variables (d-GNBP with chain variables). The d-GNBP with chain variables describes the same motion as the d-RGNBP in Equation (31), as shown by the following lemma.

Lemma 2. Suppose

• i.  
  The $(N-1)(N-2)$ vectors ${{\boldsymbol{Q}}}_{{ij}}^{(n+1)}$ and ${{\boldsymbol{P}}}_{{ij}}^{(n+1)}$ $(1\leqslant i\lt i+2\leqslant j\leqslant N)$ are given by Equations (41) and (46).
• ii.  
  The $2(N-1)$ vectors ${{\boldsymbol{Q}}}_{k,k+1}^{(n+1)}$ and ${{\boldsymbol{P}}}_{k,k+1}^{(n+1)}$ $(1\leqslant k\;\leqslant N-1)$ are the solutions of the d-GNBP with chain variables composed of Equations (47) and (48).

Then, the $N(N-1)$ vectors ${{\boldsymbol{Q}}}_{{ij}}^{(n+1)}$ and ${{\boldsymbol{P}}}_{{ij}}^{(n+1)}$ $(1\leqslant i\;\lt j\leqslant N)$ satisfy the d-RGNBP in Equation (31).

Proof. (a) Derivation of Equation (31(a)).

The vectors ${{\boldsymbol{P}}}_{{ij}}^{(n+1/2)}$ $(1\leqslant i\lt i+2\leqslant j\leqslant N)$ in Equation (45) satisfy Equation (44). Further, in combination with Equation (48), Equation (44) leads to

$$\begin{eqnarray}&&\displaystyle \frac{m}{8{m}_{i}{m}_{j}}\displaystyle \frac{{| {{\boldsymbol{Q}}}_{{ij}}^{(n+1)}| }^{2}+{| {{\boldsymbol{Q}}}_{{ij}}^{(n)}| }^{2}}{{| {{\boldsymbol{Q}}}_{{ij}}^{(n+1)}| }^{2}{| {{\boldsymbol{Q}}}_{{ij}}^{(n)}| }^{2}}{{\boldsymbol{P}}}_{{ij}}^{(n+1/2)}{\boldsymbol{K}}{({{\boldsymbol{Q}}}_{{ij}}^{(n+1/2)})}^{\top }\\ &&=\displaystyle \sum _{k=i}^{j-1}\displaystyle \frac{{{\boldsymbol{Q}}}_{k,k+1}^{(n+1)}-{{\boldsymbol{Q}}}_{k,k+1}^{(n)}}{{\rm{\Delta }}{t}^{(n,n+1)}}{\boldsymbol{K}}{({{\boldsymbol{Q}}}_{k,k+1}^{(n+1/2)})}^{\top },\\ &&1\leqslant i\lt i+2\leqslant j\leqslant N.\end{eqnarray} \tag{ 49 }$$

Similarly, the vectors ${{\boldsymbol{Q}}}_{{ij}}^{(n+1)}$ $(1\leqslant i\lt i+2\leqslant j\leqslant N)$ defined by Equation (41) fulfill Equation (40) with $l=n+1$ and ${{\boldsymbol{Q}}}_{{ij}}^{(n)}$ satisfy Equation (40) with l = n, so Equation (43) is also fulfilled. Because the right-hand side of Equation (49) coincides with that of Equation (43), Equation (31(a)) in the d-RGNBP is given.

(b) Derivation of Equation (31(b)).

Substitution of Equations (41), (45), and (46) into Equation (31(b)) yields Equation (47). Namely, Equation (47), which is part of condition (ii), equals Equation (31(b)) under condition (i). Accordingly, ${{\boldsymbol{Q}}}_{{ij}}^{(n+1)}$ and ${{\boldsymbol{P}}}_{{ij}}^{(n+1)}$ $(1\leqslant i\lt j\leqslant N)$ satisfying conditions (i) and (ii) follow from Equation (31(b)).

(c) Derivation of Equation (31(c)).

We have already stated that Equation (41) is the solution of Equation (40). Substitution of Equation (40) into the left-hand side of Equation (31c) yields $0$. Thus, ${{\boldsymbol{Q}}}_{{ij}}^{(n+1)}$ $(1\leqslant i\lt j\leqslant N)$ satisfy Equation (31(c)) under conditions (i) and (ii).□

Lemma 2 clarifies that the d-GNBP with chain variables reproduces the motion of the d-RGNBP, and Theorem 1 shows that the d-RGNBP preserves all the conserved quantities except the angular momentum. Therefore, the conservation of quantities in the d-GNBP with chain variables is stated in the following theorem.

Theorem 2. (Conserved quantities of d-GNBP with chain variables.) The d-GNBP with chain variables composed of Equations (47) and (48) exactly preserves the following three conserved quantities:

• 1.  
  the value of the Hamiltonian ${h}_{\mathrm{rel}}$ defined by Equation (9),
• 2.  
  the linear momentum ${\boldsymbol{l}}\equiv {\displaystyle \sum }_{i=1}^{N}{{\boldsymbol{p}}}_{i}^{\prime }=0$, and
• 3.  
  the position of the center of mass ${\boldsymbol{c}}\equiv {\displaystyle \sum }_{i=1}^{N}{m}_{i}{{\boldsymbol{q}}}_{i}^{\prime }=0$.

The d-GNBP with chain variables composed of Equations (47) and (48) is still singular at each ${{\boldsymbol{Q}}}_{k,k+1}^{(n)}=0$, ${{\boldsymbol{Q}}}_{k,k+1}^{(n+1)}=0$ $(1\leqslant k\leqslant N-1)$ and ${\tilde{{\boldsymbol{Q}}}}_{i,j}^{(n)}=0$, ${\tilde{{\boldsymbol{Q}}}}_{i,j}^{(n+1)}=0$ $(1\leqslant i\lt i+2\leqslant j\leqslant N)$. Therefore, the discrete-time system cannot yield regular results for collision orbits.

Using the time transformation

$$\begin{eqnarray*} & & {\rm{\Delta }}{s}^{(n,n+1)}\equiv {s}^{(n+1)}-{s}^{(n)}\\ & & =(2{\rm{\Delta }}{t}^{(n,n+1)})/\\ & & ({(\displaystyle \sum _{i=1}^{N-1}\displaystyle \frac{{m}_{i}{m}_{i+1}}{{| {{\boldsymbol{Q}}}_{i,i+1}^{(n+1)}| }^{2}}+\displaystyle \sum _{i=1}^{N-2}\displaystyle \sum _{j=i+2}^{N}\displaystyle \frac{{m}_{i}{m}_{j}}{{| {\tilde{{\boldsymbol{Q}}}}_{{ij}}^{(n+1)}| }^{2}})}^{-1}\\ & & +{(\displaystyle \sum _{i=1}^{N-1}\displaystyle \frac{{m}_{i}{m}_{i+1}}{{| {{\boldsymbol{Q}}}_{i,i+1}^{(n)}| }^{2}}+\displaystyle \sum _{i=1}^{N-2}\displaystyle \sum _{j=i+2}^{N}\displaystyle \frac{{m}_{i}{m}_{j}}{{| {\tilde{{\boldsymbol{Q}}}}_{{ij}}^{(n)}| }^{2}})}^{-1}),\end{eqnarray*}$$

we can rewrite the d-GNBP with chain variables composed of Equations (47) and (48) as the following regularized discrete-time system:

$$\begin{eqnarray}&&\displaystyle \frac{{{\boldsymbol{Q}}}_{k,k+1}^{(n+1)}-{{\boldsymbol{Q}}}_{k,k+1}^{(n)}}{{\rm{\Delta }}{s}^{(n,n+1)}}=\displaystyle \frac{m}{16{m}_{k}{m}_{k+1}}\displaystyle \frac{{| {{\boldsymbol{Q}}}_{k,k+1}^{(n+1)}| }^{2}+{| {{\boldsymbol{Q}}}_{k,k+1}^{(n)}| }^{2}}{{| {{\boldsymbol{Q}}}_{k,k+1}^{(n+1)}| }^{2}{| {{\boldsymbol{Q}}}_{k,k+1}^{(n)}| }^{2}}\\ &&\times ({(\displaystyle \sum _{i=1}^{N-1}\displaystyle \frac{{m}_{i}{m}_{i+1}}{{| {{\boldsymbol{Q}}}_{i,i+1}^{(n+1)}| }^{2}}+\displaystyle \sum _{i=1}^{N-2}\displaystyle \sum _{j=i+2}^{N}\displaystyle \frac{{m}_{i}{m}_{j}}{{| {\tilde{{\boldsymbol{Q}}}}_{{ij}}^{(n+1)}| }^{2}})}^{-1}\\ &&+{(\displaystyle \sum _{i=1}^{N-1}\displaystyle \frac{{m}_{i}{m}_{i+1}}{{| {{\boldsymbol{Q}}}_{i,i+1}^{(n)}| }^{2}}+\displaystyle \sum _{i=1}^{N-2}\displaystyle \sum _{j=i+2}^{N}\displaystyle \frac{{m}_{i}{m}_{j}}{{| {\tilde{{\boldsymbol{Q}}}}_{{ij}}^{(n)}| }^{2}})}^{-1}){{\boldsymbol{P}}}_{k,k+1}^{(n+1/2)},\\ &&1\leqslant k\leqslant N-1,\end{eqnarray} \tag{ 50a }$$

$$\begin{eqnarray}&&\displaystyle \sum _{i=1}^{j-2}{\tilde{{\boldsymbol{g}}}}_{{ij}}^{(n+1)}-\displaystyle \sum _{i=j+2}^{N}{\tilde{{\boldsymbol{g}}}}_{{ji}}^{(n+1)}=\\ &&\quad -{{\boldsymbol{g}}}_{j-1,j}^{(n+1)}+{{\boldsymbol{g}}}_{j,j+1}^{(n+1)},\quad 2\leqslant j\leqslant N-1,\end{eqnarray} \tag{ 50b }$$

$$\begin{eqnarray}&&\displaystyle \sum _{i=1}^{N-2}{\tilde{{\boldsymbol{g}}}}_{i,N}^{(n+1)}=-{{\boldsymbol{g}}}_{N-1,N}^{(n+1)},\end{eqnarray} \tag{ 50c }$$

where

$$\begin{eqnarray*}&&{{\boldsymbol{g}}}_{k,k+1}^{(n+1)}\equiv \displaystyle \frac{1}{2{| {{\boldsymbol{Q}}}_{k,k+1}^{(n+1/2)}| }^{2}}(\displaystyle \frac{{{\boldsymbol{P}}}_{k,k+1}^{(n+1)}-{{\boldsymbol{P}}}_{k,k+1}^{(n)}}{{\rm{\Delta }}{s}^{(n,n+1)}}\\ &&-\displaystyle \frac{1}{2{| {{\boldsymbol{Q}}}_{k,k+1}^{(n+1)}| }^{2}{| {{\boldsymbol{Q}}}_{k,k+1}^{(n)}| }^{2}}\\ &&\times \;({(\displaystyle \sum _{i=1}^{N-1}\displaystyle \frac{{m}_{i}{m}_{i+1}}{{| {{\boldsymbol{Q}}}_{i,i+1}^{(n+1)}| }^{2}}+\displaystyle \sum _{i=1}^{N-2}\displaystyle \sum _{j=i+2}^{N}\displaystyle \frac{{m}_{i}{m}_{j}}{{| {\tilde{{\boldsymbol{Q}}}}_{{ij}}^{(n+1)}| }^{2}})}^{-1}\\ &&+\;{(\displaystyle \sum _{i=1}^{N-1}\displaystyle \frac{{m}_{i}{m}_{i+1}}{{| {{\boldsymbol{Q}}}_{i,i+1}^{(n)}| }^{2}}+\displaystyle \sum _{i=1}^{N-2}\displaystyle \sum _{j=i+2}^{N}\displaystyle \frac{{m}_{i}{m}_{j}}{{| {\tilde{{\boldsymbol{Q}}}}_{{ij}}^{(n)}| }^{2}})}^{-1})\end{eqnarray*}$$

$$\begin{eqnarray*}&&\times (\displaystyle \frac{m}{8{m}_{k}{m}_{k+1}}({| {{\boldsymbol{P}}}_{k,k+1}^{(n+1)}| }^{2}+{| {{\boldsymbol{P}}}_{k,k+1}^{(n)}| }^{2})-2{m}_{k}{m}_{k+1}){{\boldsymbol{Q}}}_{k,k+1}^{(n+1/2)})\\ &&{\boldsymbol{K}}{({{\boldsymbol{Q}}}_{k,k+1}^{(n+1/2)})}^{\top },1\leqslant k\leqslant N-1,\\ &&{\tilde{{\boldsymbol{g}}}}_{{ij}}^{(n+1)}\equiv \displaystyle \frac{1}{2{| {\tilde{{\boldsymbol{Q}}}}_{{ij}}^{(n+1/2)}| }^{2}}(\displaystyle \frac{{\tilde{{\boldsymbol{P}}}}_{{ij}}^{(n+1)}-{\tilde{{\boldsymbol{P}}}}_{{ij}}^{(n)}}{{\rm{\Delta }}{s}^{(n,n+1)}}-\displaystyle \frac{1}{2{| {\tilde{{\boldsymbol{Q}}}}_{{ij}}^{(n+1)}| }^{2}{| {\tilde{{\boldsymbol{Q}}}}_{{ij}}^{(n)}| }^{2}}\\ &&\times \;({(\displaystyle \sum _{i=1}^{N-1}\displaystyle \frac{{m}_{i}{m}_{i+1}}{{| {{\boldsymbol{Q}}}_{i,i+1}^{(n+1)}| }^{2}}+\displaystyle \sum _{i=1}^{N-2}\displaystyle \sum _{j=i+2}^{N}\displaystyle \frac{{m}_{i}{m}_{j}}{{| {\tilde{{\boldsymbol{Q}}}}_{{ij}}^{(n+1)}| }^{2}})}^{-1}\\ &&+{(\displaystyle \sum _{i=1}^{N-1}\displaystyle \frac{{m}_{i}{m}_{i+1}}{{| {{\boldsymbol{Q}}}_{i,i+1}^{(n)}| }^{2}}+\displaystyle \sum _{i=1}^{N-2}\displaystyle \sum _{j=i+2}^{N}\displaystyle \frac{{m}_{i}{m}_{j}}{{| {\tilde{{\boldsymbol{Q}}}}_{{ij}}^{(n)}| }^{2}})}^{-1})\\ &&\times (\displaystyle \frac{m}{8{m}_{i}{m}_{j}}({| {\tilde{{\boldsymbol{P}}}}_{{ij}}^{(n+1)}| }^{2}+{| {\tilde{{\boldsymbol{P}}}}_{{ij}}^{(n)}| }^{2})-2{m}_{i}{m}_{j}){\tilde{{\boldsymbol{Q}}}}_{{ij}}^{(n+1/2)})\\ &&{\boldsymbol{K}}{({\tilde{{\boldsymbol{Q}}}}_{{ij}}^{(n+1/2)})}^{\top },1\leqslant i\lt i+2\leqslant j\leqslant N.\end{eqnarray*}$$

We call the discrete-time system in Equation (50) the discrete-time chain regularization of the GNBP (d-CRGNBP). It is clear that the d-CRGNBP in Equation (50) maintains exactly the values of the Hamiltonian, linear momentum, and position of the center of mass in Theorem 2 as does the d-GNBP with chain variables composed of Equations (47) and (48).

Because the discrete-time system in Equation (50) is regular for $| {{\boldsymbol{Q}}}_{k,k+1}^{(n)}| \to 0$ or $| {{\boldsymbol{Q}}}_{k,k+1}^{(n+1)}| \to 0,k=1,2,\cdots ,N-1$, the mass ${m}_{k}$ can have collisions or close encounters with either ${m}_{k-1}$ or ${m}_{k+1}$ without causing numerical difficulties in the time interval ${I}^{(n)}=[{t}^{(n)},{t}^{(n+1)}]$. Although the case $| {\tilde{{\boldsymbol{Q}}}}_{i,j}^{(n)}| {\text{}}{or}| {\tilde{{\boldsymbol{Q}}}}_{i,j}^{(n+1)}| (1\leqslant i\lt i+2\leqslant j\leqslant N)$ is still singular, collisions between ${m}_{i}$ and ${m}_{j}$ may be studied using the set of chained position vectors ${{\boldsymbol{Q}}}_{i,i+1}^{(n)},{{\boldsymbol{Q}}}_{i+1,i+2}^{(n)},\cdots ,{{\boldsymbol{Q}}}_{j-1,j}^{(n)}$ or the set of ${{\boldsymbol{Q}}}_{i,i+1}^{(n+1)},{{\boldsymbol{Q}}}_{i+1,i+2}^{(n+1)},\cdots ,{{\boldsymbol{Q}}}_{j-1,j}^{(n+1)}$. For $l=n,n+1,| {\tilde{{\boldsymbol{Q}}}}_{i,j}^{(l)}| \geqslant {\rm{min}}({{\boldsymbol{Q}}}_{i,i+1}^{(l)},{{\boldsymbol{Q}}}_{i+1,i+2}^{(l)},\cdots ,{{\boldsymbol{Q}}}_{j-1,j}^{(l)})$ in the chain concept of Mikkola & Aarseth (1993), so the term

$$\begin{eqnarray*}&&\displaystyle \frac{1}{2| {\tilde{{\boldsymbol{Q}}}}_{{ij}}^{(n+1)}| | {\tilde{{\boldsymbol{Q}}}}_{{ij}}^{(n)}| }({(\displaystyle \sum _{i=1}^{N-1}\displaystyle \frac{{m}_{i}{m}_{i+1}}{{| {{\boldsymbol{Q}}}_{i,i+1}^{(n+1)}| }^{2}}+\displaystyle \sum _{i=1}^{N-2}\displaystyle \sum _{j=i+2}^{N}\displaystyle \frac{{m}_{i}{m}_{j}}{{| {\tilde{{\boldsymbol{Q}}}}_{{ij}}^{(n+1)}| }^{2}})}^{-1}\\ &&+{(\displaystyle \sum _{i=1}^{N-1}\displaystyle \frac{{m}_{i}{m}_{i+1}}{{| {{\boldsymbol{Q}}}_{i,i+1}^{(n)}| }^{2}}+\displaystyle \sum _{i=1}^{N-2}\displaystyle \sum _{j=i+2}^{N}\frac{{m}_{i}{m}_{j}}{{| {\tilde{{\boldsymbol{Q}}}}_{{ij}}^{(n)}| }^{2}})}^{-1})\end{eqnarray*}$$

is numerically well behaved even for small values of $| {\tilde{{\boldsymbol{Q}}}}_{i,j}^{(n)}|$ and $| {\tilde{{\boldsymbol{Q}}}}_{i,j}^{(n+1)}|$.

In the next section, we will clarify that the d-CRGNBP can precisely reproduce some periodic orbits and non-periodic orbits that have many close encounters.

First, let us show that the the d-GNBP and d-CRGNBP methods preserve the value of the Hamiltonian ${h}_{\mathrm{rel}}$ exactly and the angular momentum ${\boldsymbol{J}}=\displaystyle \sum _{i=1}^{4}\;{\boldsymbol{p}}{^{\prime} }_{i}\times {\boldsymbol{q}}{^{\prime} }_{i}$ approximately. Figure 1 shows the dependence of the relative error growth of the values of ${h}_{\mathrm{rel}}$ and ${\boldsymbol{J}}$ for the G4BP on the above integrators. The adopted initial conditions are as follows:

$$\begin{eqnarray}&&{m}_{1}={m}_{2}={m}_{3}={m}_{4}=0.25,\\ &&{\boldsymbol{p}}{^{\prime} }_{1}=(0,0.1),\;{\boldsymbol{p}}{^{\prime} }_{2}=(0,-0.05),\\ &&{\boldsymbol{p}}{^{\prime} }_{3}=(0,-0.1),{\boldsymbol{p}}{^{\prime} }_{4}=(0,0.05),\\ &&{\boldsymbol{q}}{^{\prime} }_{1}=(-10,0),\ {\boldsymbol{q}}{^{\prime} }_{2}=(-12,0),\\ &&{\boldsymbol{q}}{^{\prime} }_{3}=(10,0),\ {\boldsymbol{q}}{^{\prime} }_{4}=(12,0).\end{eqnarray} \tag{ 51 }$$

For the above two non-regularized methods, the SI8 and d-GNBP methods, the real time step size is fixed as ${\rm{\Delta }}{t}^{(n,n+1)}=0.1$. In addition, the fictitious time step size is set to ${\rm{\Delta }}{s}^{(n,n+1)}=0.0045$, which corresponds to the average real time step size $\overline{({\rm{\Delta }}t)}=$ 0.0940412562869442696, for the leapfrog for the chain system method and ${\rm{\Delta }}{s}^{(n,n+1)}=0.0000042$, which corresponds to $\overline{({\rm{\Delta }}t)}=0.101481391749539464$, for the d-CRGNBP method. The initial condition corresponds to the Caledonian symmetric four-body problem (CS4BP), in which the configuration of the four bodies is always a parallelogram (Széll et al. 2004). Each method always sets the linear momentum ${\boldsymbol{l}}$ and the center of mass ${\boldsymbol{c}}$ at the origin in the barycentric frame.

Zoom In Zoom Out Reset image size

The relative error of the value of ${h}_{\mathrm{rel}}$ is bounded by the value ${10}^{-7}$ for the leapfrog for the chain system method, whereas the error is bounded by a sufficiently small value of ${10}^{-14}$ for the SI8, d-GNBP, and d-CRGNBP methods. In particular, the d-GNBP and d-CRGNBP methods conserve ${h}_{\mathrm{rel}}$ with an accuracy of ${10}^{-16}$. In addition, the error of ${\boldsymbol{J}}$ is bounded by 5 $\times {10}^{-5}$ for all the methods. Specifically, the SI8 and leapfrog for the chain system methods precisely preserve the value of ${\boldsymbol{J}}$ with an accuracy of ${10}^{-16}$. Because no integration method preserves both ${h}_{\mathrm{rel}}$ and ${\boldsymbol{J}}$, none would seem to reproduce stable orbits for a long time interval. However, Section 4.2 will clarify that this prediction is incorrect.

In Section 4.1, we showed that the d-GNBP and d-CRGNBP methods do not accurately preserve the angular momentum. Therefore, it is questionable whether they can reproduce various orbits of the GNBP because the orbits lie on the manifold determined by the conserved quantities. To answer this question, we show that some periodic orbits computed by these methods accurately coincide with those of the G3BP, G4BP, and G5BP.

4.2.1. Periodic Orbits in the G4BP

First, we show that the d-CRGNBP can precisely follow the orbits for a family of stable periodic orbits in the G4BP. Yan & Ouyang (2014) introduced a variational method to provide some stable periodic orbits in the G3BP, G4BP, and G5BP, and then they listed their initial conditions. We rewrite one of the initial conditions for stable periodic orbits in the G4BP, which corresponds to Figure 25 of Yan & Ouyang (2014), as follows:

$$\begin{eqnarray}&&{m}_{1}={m}_{2}={m}_{3}={m}_{4}=1,\\ &&{\boldsymbol{p}}{^{\prime} }_{1}=(-0.6978,-0.5025),\ {\boldsymbol{p}}{^{\prime} }_{2}=(-0.6978,0.5025),\\ &&{\boldsymbol{p}}{^{\prime} }_{3}=(0.0178,0),{\boldsymbol{p}}{^{\prime} }_{4}=(1.3777,0),\\ &&{\boldsymbol{q}}{^{\prime} }_{1}=(0.6675,1.1148),\ {\boldsymbol{q}}{^{\prime} }_{2}=(-0.6675,1.1148),\\ &&{\boldsymbol{q}}{^{\prime} }_{3}=(0,-1.7197),\ {\boldsymbol{q}}{^{\prime} }_{4}=(0,-0.5099).\end{eqnarray} \tag{ 52 }$$

We introduce a uniform rotating frame $O-{x}_{[1]}{x}_{[2]}$ in which the origin stays at the center of mass, and the ${x}_{[1]}$ axis passes through the origin and the mass ${m}_{1}$. We integrated over the time interval 0 $\leqslant t\leqslant 10,000$ and used the common fixed real time step ${\rm{\Delta }}{t}^{(n,n+1)}=0.1$ for the SI8 and d-GNBP methods, the fixed fictitious time step ${\rm{\Delta }}{s}^{(n,n+1)}=0.33$, which corresponds to the average time step size $\bar{({\rm{\Delta }}t)}=0.0963516604110964389$, for the leapfrog for the chain system method, and ${\rm{\Delta }}{s}^{(n,n+1)}=0.017$, which corresponds to $\bar{({\rm{\Delta }}t)}=0.101449146855772802$, for the d-CRGNBP method. Each orbit computed by the leapfrog for chain system and d-CRGNBP methods stays near a closed periodic orbit for a long time, whereas the SI8 and d-GNBP methods yield no closed orbit. Figure 2 shows the orbits of all four particles in the $O-{x}_{[1]}{x}_{[2]}$ frame for only the SI8 and d-CRGNBP methods. Note that the leapfrog for the chain system method yields four orbits, each of which is similar to the orbit computed by the d-CRGNBP method, and that the d-GNBP and SI8 methods yield no closed orbit. The result would mean that only regularized integration methods can reproduce these stable periodic orbits for a long time interval.

Zoom In Zoom Out Reset image size

4.2.2. Periodic Orbits in the G5BP

The restricted four-body problem (R4BP), like the R3BP, has some periodic orbits, the orbits of elliptic relative equilibrium solutions (Majorana 1981; Baltagiannis & Papadakis 2011a; Minesaki 2015). In these orbits, each of the three primaries, ${m}_{1},{m}_{2}$, and ${m}_{3}$, is fixed at the apex of an equilateral triangle that rotates while alternately expanding and contracting. Namely, these primaries revolve in the orbits of an elliptic Lagrangian triangle solution. The massless particle ${m}_{4}$ does not influence the primaries ${m}_{1},{m}_{2}$, and ${m}_{3}$. Further, ${m}_{4}$ and these primaries form a tetragon whose shape is invariant. Each orbit of ${m}_{1},{m}_{2},{m}_{3}$, and ${m}_{4}$ is a conic curve in the inertial barycentric frame.

In the $2+3$-body problem, two massless particles, ${m}_{4}$ and ${m}_{5}$, move in the gravitational field of three primaries, ${m}_{1},{m}_{2}$, and ${m}_{3}$, in the periodic orbits of an elliptic Lagrangian triangle solution. In this problem, each of the five particles, ${m}_{1},{m}_{2},{m}_{3},{m}_{4}$, and ${m}_{5}$, is fixed at a pentagon whose shape is invariant. Because one massless particle, ${m}_{4}$ or ${m}_{5}$, does not influence the three primaries ${m}_{1},{m}_{2}$, and ${m}_{3}$ or the other massless particle, each orbit of ${m}_{4}$ and ${m}_{5}$ coincides with one of two different orbits of ${m}_{4}$ in the R4BP.

In this section, we clarify that the d-CRGNBP accurately computes the orbits of some stable periodic solutions in the G5BP that reduce to the above orbits in the $2+3$-body problem as the masses ${m}_{4}$ and ${m}_{5}$ go to zero. The initial conditions are given for one such periodic solution in the G5BP as

$$\begin{eqnarray}&&{m}_{1}=0.9981455689622930,\\ &&{m}_{2}={m}_{3}=0.0009272155188535,\\ &&{m}_{4}={m}_{5}=1.0\times {10}^{-14},\\ &&{{\boldsymbol{p}}}_{1}^{\prime }=(0.0003993831619767,0.0015524567530639),\\ &&{{\boldsymbol{p}}}_{2}^{\prime }=(-0.0006486798609336,-0.0006607221902365),\\ &&{{\boldsymbol{p}}}_{3}^{\prime }=(0.0002492966989569,-0.0008917345628274),\\ &&{{\boldsymbol{p}}}_{4}^{\prime }=(-9.991644458074203\\ &&\qquad \times {10}^{-15},2.30012936560301\times {10}^{-16}),\\ &&{{\boldsymbol{p}}}_{5}^{\prime }=(8.862203117935508\times {10}^{-15},\\ &&\qquad -4.620305273897407\times {10}^{-15}),\\ &&{{\boldsymbol{q}}}_{1}^{\prime }=(0.0013908232782802,0.0008029921941103),\\ &&{{\boldsymbol{q}}}_{2}^{\prime }=(-0.9986091767217197,0.0008029921941103),\\ &&{{\boldsymbol{q}}}_{3}^{\prime }=(-0.4986091767217198,-0.8652224115903283),\\ &&{{\boldsymbol{q}}}_{4}^{\prime }=(-0.6829889028243268,0.7296470433670392),\\ &&{{\boldsymbol{q}}}_{5}^{\prime }=(0.2903984239398986,-0.9563092620322480).\end{eqnarray} \tag{ 53 }$$

As ${m}_{4}\to$ 0, the values of ${\boldsymbol{p}}{^{\prime} }_{4}$ and ${\boldsymbol{q}}{^{\prime} }_{4}$ reduce to the initial conditions for a stable relative equilibrium orbit corresponding to one equilibrium point ${L}_{5}$ in the R4BP (see the values of ${\boldsymbol{q}}{^{\prime} }_{4}$ and ${\boldsymbol{v}}{^{\prime} }_{4}$ in Equation (45) of Minesaki 2015). In addition, in the limit of ${m}_{5}\to$ 0, the values of ${\boldsymbol{p}}{^{\prime} }_{5}$ and ${\boldsymbol{q}}{^{\prime} }_{5}$ lead to the initial conditions for a stable relative equilibrium orbit corresponding to the other equilibrium point ${L}_{6}$ in the R4BP (see the values of ${\boldsymbol{q}}{^{\prime} }_{4}$ and ${\boldsymbol{v}}{^{\prime} }_{4}$ in Equation (46) of Minesaki 2015).

The common real time steps for the SI8 and d-GNBP methods were fixed at ${\rm{\Delta }}{t}^{(n,n+1)}=0.1$, and we set the fixed fictitious time steps ${\rm{\Delta }}{s}^{(n,n+1)}=0.00019$, which corresponds to the average time step size $\overline{({\rm{\Delta }}t)}=0.102601242925268408$, for the leapfrog for the chain system method and ${\rm{\Delta }}{s}^{(n,n+1)}=0.058$, which corresponds to $\overline{({\rm{\Delta }}t)}=0.103974428092040402$, for the d-CRGNBP. We integrated over the time interval 0 $\leqslant \;t\;\leqslant$ 1000. The d-GNBP or d-CRGNBP method reproduces closed periodic orbits over long time intervals, whereas the SI8 and the leapfrog for the chain system methods do not yield even the three elliptic orbits of the primaries ${m}_{1},{m}_{2}$, and ${m}_{3}$. Figure 3 shows the orbits of all five particles in the barycentric frame for only the SI8 and d-GNBP methods. The result suggests that only energy–momentum methods such as the d-GNBP and d-CRGNBP methods can precisely trace these stable periodic orbits for a long time interval.

Zoom In Zoom Out Reset image size

4.2.3. Periodic Orbits in the G3BP

As described in Sections 4.2.1 and 4.2.2, only the d-CRGNBP method can follow both periodic orbits in the G4BP and G5BP. In this section, we now clarify that the d-CRGNBP method can also precisely reproduce some stable periodic orbits in the G3BP. Tsouroplis & Zagouras (1984) numerically provided various orbits of asymmetric periodic solutions for the planar G3BP and confirmed their stability. The adopted initial conditions are the first initial conditions in Table V of Tsouroplis & Zagouras (1984), as follows:

$$\begin{eqnarray}&&{m}_{1}={m}_{2}=0.4974,\ {m}_{3}=0.0052,\\ &&{{\boldsymbol{p}}}_{1}^{\prime }=(0.0002239014509303,-0.2597669231723695),\\ &&{{\boldsymbol{p}}}_{2}^{\prime }=(-0.0004670997030698,0.2446876785996305),\\ &&{{\boldsymbol{p}}}_{3}^{\prime }=(0.0002431982521395,0.0150792445727390),\\ &&{{\boldsymbol{q}}}_{1}^{\prime }=(-0.506855165842,0),\ {\boldsymbol{q}}{^{\prime} }_{2}=(0.512272834158,0),\\ &&{{\boldsymbol{q}}}_{3}^{\prime }=(-0.267445165842,0).\end{eqnarray} \tag{ 54 }$$

These conditions correspond to a set of stable periodic orbits.

We apply the SI8, GNBP, leapfrog for the chain system, and d-CRGNBP methods to the G3BP with the initial conditions in Equation (54) and integrate over the time interval 0 $\leqslant \;t\;\leqslant$ 1000. The real time step size is set to ${\rm{\Delta }}{t}^{(n,n+1)}=0.01$ for the SI8 and GNBP methods, and the fictitious time step sizes are fixed as ${\rm{\Delta }}{s}^{(n,n+1)}=0.002821$, which corresponds to the average time step size $\overline{({\rm{\Delta }}t)}=0.008230723645$, for the leapfrog for the chain system method and ${\rm{\Delta }}{s}^{(n,n+1)}=0.04$, which corresponds to $\overline{({\rm{\Delta }}t)}=0.009572767391$, for the d-CRGNBP. Figure 4 shows the orbits of ${m}_{1},{m}_{2}$, and ${m}_{3}$ obtained by the leapfrog for the chain system and d-CRGNBP methods in the coordinate system $O\;-\;{x}_{[1]}{x}_{[2]}$ described in Section 4.2.1.

Zoom In Zoom Out Reset image size

Theoretically, the mass ${m}_{1}$ moves on the ${x}_{[1]}$ axis, and the masses ${m}_{2}$ and ${m}_{3}$ each revolve in a closed orbit. For the time interval 0 $\leqslant \;t\leqslant 10,000$, the SI8, leapfrog for the chain system, and d-GNBP methods do not yield closed orbits, and only the d-CRGNBP method accurately reproduces them. The results obtained by the SI8 and leapfrog for the chain system methods show that each particle escapes after staying near a closed orbit for a few periods, and those given by the d-GNBP method indicate that all the particles escape without staying near closed orbits. (However, Figure 4 does not include the results for the SI8 and d-GNBP methods.) This means that only the d-CRGNBP can accurately compute the periodic orbits in the G3BP.

We clarify that the d-CRGNBP method computes various stable periodic orbits in the G3BP, G4BP, and G5BP more accurately than the geometric integration methods (the SI8, d-GNBP, and leapfrog for the chain system methods).

To avoid significant rounding errors caused by the singularities in the GNBP, the chain concept of Mikkola & Aarseth (1993) and adequate time transformations are introduced for the leapfrog for the chain system and d-CRGNBP methods. To confirm their efficacy, the leapfrog for the chain system and d-CRGNBP methods are integrated through many close encounters from $t=0$ to about a given real time $t=T$ and then backward to about $t=0$. For each particle ${m}_{i}$, we check whether its computed orbit from $t=0$ to about $t=T$ coincides with its integrated orbit backward to about $t=0$. We call this approach a reversibility check.

We give the initial conditions for a test G5BP:

$$\begin{eqnarray}&&{m}_{1}=0.992,{m}_{2}={m}_{3}={m}_{4}={m}_{5}=0.002,\\ &&{{\boldsymbol{p}}}_{1}^{\prime }=(0,0),\ {\boldsymbol{p}}{^{\prime} }_{2}=(0,0.02),{\boldsymbol{p}}{^{\prime} }_{3}=(0,-0.01),\\ &&{{\boldsymbol{p}}}_{4}^{\prime }=(0,-0.02),\ {\boldsymbol{p}}{^{\prime} }_{5}=(0,0.01),\\ &&{{\boldsymbol{q}}}_{1}^{\prime }=(0,0),\ {\boldsymbol{q}}{^{\prime} }_{2}=(-0.015,0),\ {\boldsymbol{q}}{^{\prime} }_{3}=(-0.03,0),\\ &&{{\boldsymbol{q}}}_{4}^{\prime }=(0.015,0),\ {\boldsymbol{q}}{^{\prime} }_{5}=(0.03,0),\end{eqnarray} \tag{ 55 }$$

and set T = 0.7. Theoretically, these conditions correspond to the Caledonian symmetric five-body problem (CS5BP) in which one body is fixed at the origin, and the other four bodies are always configured in a parallelogram (Shoaib et al. 2013).

We define the functions ${r}_{\mathrm{min}}(t)$ and ${r}_{\mathrm{max}}(t)$ as the minimum and maximum of $| {{\boldsymbol{q}}}_{{ij}}| (1\leqslant i\lt j\leqslant 5)$ at a real time $t$, respectively. In this test, two-body close encounters occur 12 times until $t=T=0.7$, as shown in Figure 5. We consider that close encounters occur when ${r}_{\mathrm{min}}(t)/{r}_{\mathrm{max}}(t)\lt 0.05$.

Zoom In Zoom Out Reset image size

We show the results of this check for the leapfrog for the chain system and d-CRGNBP methods. We used the fictitious time step size ${\rm{\Delta }}{s}^{(n,n+1)}=0.0000099$, which corresponds to the average real time step size ${\rm{\Delta }}{t}^{(n,n+1)}=0.0000335496683762689432$, for the leapfrog for the chain system method and ${\rm{\Delta }}{s}^{(n,n+1)}=0.00001$, which corresponds to ${\rm{\Delta }}{t}^{(n,n+1)}\;=0.0000338629446219376209$, for the d-CRGNBP method. Backward integration of the leapfrog method yields an orbit for each particle that differs from that computed by forward integration. In contrast, backward integration of the d-CRGNBP method traces the same orbit of each particle as forward integration. In Figure 6, the orbit of only particle ${m}_{3}$ is drawn for the leapfrog and d-CRGNBP methods.

Zoom In Zoom Out Reset image size