N-body integration

The Hermite 4th order integrator (H4) [MakinoAarseth1992] is a widely used N- body integrator scheme, which is used in the famous Aarseth integrator family [Aarseth1999] [Aarseth2003]. This integrator are useful as an starting point because:

  • First of all, all Aarseth codes are published, so it is easy to compare new implementations, with an existing working code.
  • The NBODY families are an evolution of a N- body integrator, so it is convenient to take a look of previous version, to understand determined process.
  • The publications around Aarseth codes is enough to build up an idea of how its integrator works.

The H4 integrator scheme is based on a predictor-corrector scenario, this means, that we will use an extrapolation of the equations of motion to get a predicted position and velocity at some time, then we will use this information to get the new accelerations, then corrects the predicted values using interpolation, based on finite differences terms. We can use polynomial adjustment in the gravitational forces evolution among the time, because the force acting over each particle changes smoothly.

To predict the equation of motion, we have an easy to evaluate term, which is an explicit polynomial. This prediction is less accurate, but it is improved in the corrector phase, which consist of an implicit polynomial that will require good initial values to scale to a good convergence.

This algorithm it is called fourth-order, because the predictor consider the contributions of the third order polynomial, then, after obtaining the accelerations, adds a fourth-order corrector term.

Hermite 4th order

The mathematical formulation of the integrator main steps are the following:

Prediction,

\[\begin{split}\boldsymbol{r}_{i,\rm pred} &= \boldsymbol{r}_{i,0} + \boldsymbol{v}_{i,0} \Delta t_{i} + \frac{1}{2!} \boldsymbol{a}_{i,0} \Delta t^{2}_{i} + \frac{1}{3!} \boldsymbol{\dot{a}}_{i,0} \Delta t^{3}_{i} \\ \boldsymbol{v}_{i,\rm pred} &= \boldsymbol{v}_{i,0} + \boldsymbol{a}_{i,0} \Delta t_{i} + \frac{1}{2!} \boldsymbol{\dot{a}}_{i,0} \Delta t^{2}_{i}\end{split}\]

Force calculation,

\[\begin{split}\boldsymbol{a}_{i,1} &= \sum_{\substack{j=0\\j\neq i}}^{N} G m_{j} \frac{\boldsymbol{r}_{ij}} {(r_{ij}^2 + \epsilon^{2})^{\frac{3}{2}}}, \\ \boldsymbol{\dot{a}}_{i,1} &= \sum_{\substack{j=0\\j\neq i}}^{N} G m_{j} \left[ \frac{\boldsymbol{v}_{ij}} {(r_{ij}^2 + \epsilon^{2})^{\frac{3}{2}}} - \frac{3(\boldsymbol{v}_{ij}\cdot \boldsymbol{r}_{ij}) \boldsymbol{r}_{i}} {(r_{ij}^2 + \epsilon^{2})^{\frac{5}{2}}} \right],\end{split}\]

where,

\[\begin{split}\boldsymbol{r}_{ij} &= \boldsymbol{r}_{j,pred} - \boldsymbol{r}_{i,pred}, \\ \boldsymbol{v}_{ij} &= \boldsymbol{v}_{j,pred} - \boldsymbol{v}_{i,pred}, \\ r_{ij} &= |\boldsymbol{r}_{ij}|\end{split}\]

is important to note, that \((\boldsymbol{v}_{ij}\cdot \boldsymbol{r}_{ij})\) correspond to the dot product, and not a simple multiplication.

Correction,

\[\begin{split}\boldsymbol{r}_{i,1} &= \boldsymbol{r}_{i,pred} + \frac{1}{24} \Delta t_{i}^{4} \boldsymbol{a}_{i,0}^{(2)} + \frac{1}{120} \Delta t_{i}^{5} \boldsymbol{a}_{i,0}^{(3)} \\ \boldsymbol{v}_{i,1} &= \boldsymbol{v}_{i,pred} + \frac{1}{4} \Delta t_{i}^{3} \boldsymbol{a}_{i,0}^{(2)} + \frac{1}{24} \Delta t_{i}^{4} \boldsymbol{a}_{i,0}^{(3)}\end{split}\]

Calculate the 2nd and the 3rd derivative of the acceleration \((\boldsymbol{a}_{i,1}^{(2)}, \boldsymbol{a}_{i,1}^{(3)})\) using the third-order Hermite interpolation polynomial constructed using \(\boldsymbol{a}_{i,0}\) and \(\boldsymbol{\dot{a}}_{i,0}\):

\[\begin{split}\boldsymbol{a}_{i,1}(t) &= \boldsymbol{a}_{i,0} + \boldsymbol{\dot{a}}_{i,0} \Delta{t}_{i,0} + \frac{1}{2} \Delta t^{2}_{i,0} \boldsymbol{a}_{i,0}^{(2)} + \frac{1}{6} \Delta t^{3}_{i,0} \boldsymbol{a}_{i,0}^{(3)}\\\end{split}\]

where \(\boldsymbol{a}_{i,0}\) and \(\boldsymbol{\dot{a}}_{i,0}\) are the acceleration and jerk calculated at the previous time \(t\), the second and third acceleration derivatives \(\boldsymbol{a}_{i}^{(2)}\) and \(\boldsymbol{a}_{i}^{(3)}\) are given by:

\[\begin{split}\boldsymbol{a}_{i,0}^{(2)} &= \frac{ -6 (\boldsymbol{a}_{i,0} - \boldsymbol{a}_{i,1}) - \Delta t_{i,0} (4 \boldsymbol{\dot{a}}_{i,0} + 2 \boldsymbol{\dot{a}}_{i,1})} {\Delta t_{i,0}^{2}} \\ \boldsymbol{a}_{i,0}^{(3)} &= \frac{-12 (\boldsymbol{a}_{i,0} - \boldsymbol{a}_{i,1}) - 6 \Delta t_{i,0} (\boldsymbol{\dot{a}}_{i,0} + \boldsymbol{\dot{a}}_{i,1})} {\Delta t_{i,0}^{3}}\end{split}\]

where \(\boldsymbol{a}_{i,1}\) and \(\boldsymbol{\dot{a}}_{i,1}\) are the acceleration and the jerk at the time \(t_{i} + \Delta t_{i}\).

Finally, it is necessary to calculate the next time-step for the \(i\) particle (\(\Delta t_{i,1}\)) and time \(t\) using the following formulas:

(1)\[\begin{split}t_{i,1} &= t_{i,0} + \Delta t_{i,0} \\ \Delta t_{i,1} &= \sqrt{\eta \frac{ |\boldsymbol{a}_{i,1}| |\boldsymbol{a}_{i,1}^{(2)}| + |\boldsymbol{j}_{i,1}|^{2} } {|\boldsymbol{j}_{i,1}| |\boldsymbol{a}_{i,1}^{(3)}| + |\boldsymbol{a}_{i,1}^{(2)}|^{2}}}\end{split}\]

where \(\eta\) the accuracy control parameter, \(\boldsymbol{a}_{i,1}\) and \(\boldsymbol{\dot{a}}_{i,1}\) are already known, \(\boldsymbol{a}_{i,1}^{(3)}\) has the same value as \(\boldsymbol{a}_{i,0}^{(3)}\), due the third-order interpolation; and \(\boldsymbol{a}_{i,1}^{(2)}\) is given by:

\[\begin{split}\boldsymbol{a}_{i,1}^{(2)} &= \boldsymbol{a}_{i,0}^{(2)} + \Delta t_{i}\ \boldsymbol{a}_{i,0}^{(3)}\\\end{split}\]

Block time steps

The main idea to use block time steps is to have several blocks (group) of particles sharing the same time steps, this will decrease the amount of operations of the integration process obtaining a similar accuracy than the individual time steps scheme.

In this scenario, the parallelization fits very good because compared with the individual scheme, we can have several chunks of threads working with different particles blocks.

The particle \(i\), will be part of the lower \(n-\) th block time step, between two powers of two quantities,

\[2^{n} \Delta t_{s} \leq \Delta t_{i} < 2^{n+1} \Delta t_{s}.\]

where \(\Delta t_{i}\) is determined by equation (1), and \(\Delta t_{s}\) is a constant.

The distribution of the particles among the blocks is determined by the following condition,

\[\Delta t_{i, new} = 2^{\lceil{\log_{2}{\Delta t_{i}}}\rceil - 1},\]

noindent and is described in the Figure~ref{fig:block_time steps}.

Block time-steps diagram

Block time steps illustration.

The different blocks are represented by different colours. Each particle is predicted (not move) at every time \(t\) (gray arrows), even if it’s not their block time-step (blue circles). The particles will be updated (moved) only in their block time-step (black circles). At every time that a particle is updated, they can change its block time-step, in this case, the particles \(1, 4\) and \(N\) change their block.}

For the boundaries of the time steps, we use the Aarseth criterion [Aarseth2003] for the lower limit,

\[\Delta t_{min} = 0.04 \left(\frac{\eta_{I}}{0.2}\right)^{1/2} \left(\frac{R^{3}_{\rm cl}}{\bar{m}}\right)^{1/2},\]

where \(\eta_{I}\) is the initial parameter for accuracy, which is typically \(0.01\); \(R^{3}_{\rm cl}\) is the close encounter distance; and \(\bar{m}\) is the mean mass.

Typically \(\Delta t_{\rm min} = 2^{-23}\). On the other hand, we set a maximum time step \(\Delta t_{\rm max} = 2^{-3}\). When updating a particle’s time-step, if it is out the this boundaries, we modify the value to \(\Delta t_{\rm min}\) if \(\Delta t < \Delta t_{\rm min}\), and to \(\Delta t_{\rm max}\) if \(\Delta t > \Delta t_{\rm max}\).

[MakinoAarseth1992]On a Hermite integrator with Ahmad-Cohen scheme for gravitational many-body problems http://adsabs.harvard.edu/abs/1992PASJ...44..141M
[Aarseth1999]From NBODY1 to NBODY6: The Growth of an Industry http://adsabs.harvard.edu/abs/1999PASP..111.1333A
[Aarseth2003](1, 2) Gravitational N-Body Simulations http://adsabs.harvard.edu/abs/2003gnbs.book.....A