---
author:
- Walter Dehnen
title: 'A Hierarchical ${\cal O}(N)$ Force Calculation Algorithm'
---

Introduction {#sec:intro}
============

In $N$-body simulations of stellar dynamics (or any other dynamics incorporating long-range forces), the computation, at every time step, of the gravitational forces between $N$ mutually interacting bodies dominates the operational effort. In many situations, such as studies of collisionless stellar systems, the error of these simulations is dominated by the noise in the distribution of the bodies, whose number $N$ is just a numerical parameter. Therefore, instead of computing the forces exactly by direct summation over all pairs of bodies, one may use approximate but much faster methods, allowing substantially larger $N$ and hence significantly reduced noise.

In stellar dynamics, one of the most commonly used approximate methods is the Barnes & Hut tree code [@BH] and its clones. With these methods, one first sorts the bodies into a hierarchical tree of cubic cells and pre-computes multipole moments of each cell. Next, the force at any body’s position and generated by the contents of some cell is calculated by a multipole expansion if the cell is well-separated from the body; otherwise the forces generated by the cell’s child nodes are taken. This technique reduces the number of interactions per body to ${\cal O}(\log N)$ and hence results in an overall complexity of ${\cal O}(N\log N)$.

Another technique used frequently, for instance in molecular dynamics, is Greengard & Rokhlin’s [@FMMa; @FMMa1] fast multipole method (FMM) and its variants. Traditionally, these methods first sort the bodies into a hierarchy of nested grids, pre-compute multipole moments of each cell, and then compute the forces between grid cells by a multipole expansion, usually in spherical harmonics. That is, cells are not only sources but also sinks, which formally reduces the complexity to ${\cal O}(N)$, although the author is not aware of an empirical demonstration in three dimensions[^1].

Currently, no useful implementation of FMM for application in stellar dynamics exists. In fact, it has been shown [@CM] that, for this purpose, the FMM in its traditional form cannot compete with the tree code. The reason, presumably, is two-fold: first, because stellar systems are very inhomogeneous, non-adaptive methods are less useful; second FMM codes are traditionally designed for high accuracy, whereas in collisionless stellar dynamics a relative force error of few $10^{-3}$ is often sufficient.

Here, we describe in detail a new code, designed for application in the low-accuracy regime, that combines the tree code and FMM whereby taking the better of each. In order to be fully adaptive, we use a hierarchical tree of cubic cells. The force is calculated employing [*mutual*]{} cell-cell interactions, in which both cells are source and sink simultaneously. Whether a given cell-cell interaction can be executed or must be split, is decided using an improved multipole-acceptance criterion (MAC). The new code is a further development of the code presented in [@Da], which in turn may be considered an improvement of a code given in [@WS]. It is substantially faster than the tree code and empirically shows a complexity of ${\cal O}(N)$ or even less.

The paper is organized as follows. In [§\[sec:grav\]]{}, the numerical concepts are presented; [§\[sec:alg\]]{} describes the algorithm; in [§\[sec:assess\]]{} the force errors are empirically assessed for some typical stellar dynamical test cases; [§\[sec:perform\]]{} presents empirical measurements of the complexity and performance, also in comparison to other methods; finally, [§\[sec:summ\]]{} sums up and concludes.

Approximating Gravity {#sec:grav}
=====================

The goal is to approximately compute the gravitational potential $\Phi$ and its derivative, the acceleration, at all body positions ${{\mbox{\boldmath$x$}}}_i$ and generated by all other $N-1$ bodies $$\label{potential}
\Phi({{\mbox{\boldmath$x$}}}_i) =-\sum_{j\neq
  i}\mu_j\,g({{\mbox{\boldmath$x$}}}_i-{{\mbox{\boldmath$x$}}}_j),$$ where $\mu_i$ is the weight of the $i$th body.

Consider two cells A and B, see Figure \[fig:cells\], with centers of mass ${{\mbox{\boldmath$z$}}}_{{\rm A}}$ and ${{\mbox{\boldmath$z$}}}_{{\rm B}}$, respectively. The Greens function describing the mutual interaction between a body at position ${{\mbox{\boldmath$x$}}}$ in cell A and a body at position ${{\mbox{\boldmath$y$}}}$ in cell B may be Taylor expanded about the separation ${{\mbox{\boldmath$R$}}}\equiv{{\mbox{\boldmath$z$}}}_{{\rm A}}-{{\mbox{\boldmath$z$}}}_{{\rm B}}$: $$\label{Taylor}
  g({{\mbox{\boldmath$x$}}}-{{\mbox{\boldmath$y$}}}) = \sum_{n=0}^p {1\over n!}\,({{\mbox{\boldmath$x$}}}-{{\mbox{\boldmath$y$}}}-{{\mbox{\boldmath$R$}}})^{(n)}
  \odot{\mbox{\boldmath$\nabla$}}^{(n)}g({{\mbox{\boldmath$R$}}})\;+\;{\cal R}_p(g),$$ where $p$ is the order of the expansion, while ${\cal R}_p$ denotes the Taylor series remainder (see also Appendix A). Here, we follow Warren & Salmon [@WS] by using the notational shorthand in which ${{\mbox{\boldmath$x$}}}^{(n)}$ indicates the $n$-fold outer product of the vector ${{\mbox{\boldmath$x$}}}$ with itself, while $\odot$ denotes a tensor inner product. When inserting [(\[Taylor\])]{} into [(\[potential\])]{}, whereby restricting the sum over $j$ to all bodies within cell B, we obtain for the potential at every position ${{\mbox{\boldmath$x$}}}$ in cell A and generated by all bodies in cell B [@WS] $$\begin{aligned}
 \label{local-expn}
  \Phi_{{{\rm B}}\to{{\rm A}}}({{\mbox{\boldmath$x$}}}) &=&-\sum_{m=0}^p {1\over m!}\,({{\mbox{\boldmath$x$}}}-{{\mbox{\boldmath$z$}}}_{{\rm A}})^{(m)} 
                            \odot{{\mbox{\boldmath$\sf C$}}^{m,p}}_{{{\rm B}}\to{{\rm A}}}
                            \;+\;{\cal R}_p(\Phi_{{{\rm B}}\to{{\rm A}}})
                            \\[1ex] \label{coeffs}
  {{\mbox{\boldmath$\sf C$}}^{m,p}}_{{{\rm B}}\to{{\rm A}}}  &=& \sum_{n=0}^{p-m} {(-1)^n\over n!}\;{\mbox{\boldmath$\nabla$}}^{(n+m)}
                            g({{\mbox{\boldmath$R$}}}) \odot {{\mbox{\boldmath$\sf M$}}^{n}}_{{\rm B}},
                            \\[1ex] \label{multipoles}
  {{\mbox{\boldmath$\sf M$}}^{n}}_{{\rm B}}&=& \sum_{{\mbox{\boldmath$\scriptstyley$}}_i\in{{\rm B}}} \mu_i\,({{\mbox{\boldmath$y$}}}_i-{{\mbox{\boldmath$z$}}}_{{\rm B}})^{(n)}.\end{aligned}$$ The summation over $m$ in [equation [(\[local-expn\])]{}]{} represents the evaluation of gravity, represented by the field tensors ${{\mbox{\boldmath$\sf C$}}^{m,p}}_{{{\rm B}}\to{{\rm A}}}$, at the evaluation point ${{\mbox{\boldmath$x$}}}$ within the sink cell A. The computation of the field tensors via the summation over $n$ in [equation [(\[coeffs\])]{}]{} represents the interaction between sink cell A and source cell B, represented by its multipole moments ${{\mbox{\boldmath$\sf M$}}^{n}}_{{\rm B}}$.

The symmetry between  and  of the Taylor expansion [(\[Taylor\])]{} has the important consequence that, if [equations [(\[local-expn\])]{}]{} and [(\[coeffs\])]{} are used to compute ${\mbox{\boldmath$\nabla$}}\Phi_{{{\rm B}}\to{{\rm A}}}$ [*and*]{} ${\mbox{\boldmath$\nabla$}}\Phi_{{{\rm A}}\to{{\rm B}}}$, Newton’s third law is satisfied by construction. For instance, the sum over all forces of $N$ bodies vanishes within floating point accuracy.

Note that the highest-order multipole moments, ${{\mbox{\boldmath$\sf M$}}^{n=p}}$, contribute only to the coefficients ${{\mbox{\boldmath$\sf C$}}^{0,p}}$, and hence affect only $\Phi$ but not ${\mbox{\boldmath$\nabla$}}\Phi$. Since, in stellar dynamics the acceleration rather than the potential is to be computed, one may well ignore ${{\mbox{\boldmath$\sf M$}}^{p}}$, reducing CPU-time and memory requirements.

The formulae used in the standard tree code can be obtained by setting ${{\mbox{\boldmath$x$}}}={{\mbox{\boldmath$z$}}}_{{\rm A}}$, corresponding to body sinks. In this case, potential and acceleration are approximated by $-{{\mbox{\boldmath$\sf C$}}^{0,p}}_{{{\rm B}}\to{{\rm A}}}$ and ${{\mbox{\boldmath$\sf C$}}^{1,p}}_{{{\rm B}}\to{{\rm A}}}$, respectively.

Gravity Between Well-Separated Cells {#sec:hard}
------------------------------------

In our implementation, we stick to a third order expansion ($p=3)$, whereby ignoring octopoles ${{\mbox{\boldmath$\sf M$}}^{3}}_{{\rm B}}$ (see remark above). The dipole ${{\mbox{\boldmath$\sf M$}}^{1}}_{{\rm B}}$ vanishes by construction and the Taylor-series coefficients for spherical Greens functions read (with Einstein’s sum convention) $$\label{hard-coeffs}
  \begin{array}{lcl}
    {{\sf C}^{0}}_{{{\rm B}}\to{{\rm A}}}        &=& {{\sf M}^{}}_{{\rm B}}\big[D^0
                                  + {{\textstyle{1\over2}}}{\tilde{\sf M}{}^{2}}_{{{\rm B}}ii}\,D^1
                                  + {{\textstyle{1\over2}}}R_iR_j{\tilde{\sf M}{}^{2}}_{{{\rm B}}ij}D^2
                                    \big],                           \\[1ex] 
    {{\sf C}^{1}}_{{{\rm B}}\to{{\rm A}},\,i}    &=& {{\sf M}^{}}_{{\rm B}}\big[R_i\big(D^1      
                                  + {{\textstyle{1\over2}}}{\tilde{\sf M}{}^{2}}_{{{\rm B}}jj}D^2
                                  + {{\textstyle{1\over2}}}R_jR_k{\tilde{\sf M}{}^{2}}_{{{\rm B}}jk}D^3\big)
                                  + R_j{\tilde{\sf M}{}^{2}}_{{{\rm B}}ij}D^2
                                    \big],                           \\[1ex]
    {{\sf C}^{2}}_{{{\rm B}}\to{{\rm A}},\,ij}   &=& {{\sf M}^{}}_{{\rm B}}\big[
                                    \delta_{ij}\,D^1+ R_iR_j\,D^2
                                    \big],                           \\[1ex]
    {{\sf C}^{3}}_{{{\rm B}}\to{{\rm A}},\,ijk}  &=& {{\sf M}^{}}_{{\rm B}}\big[\big(
                                    \delta_{ij}R_k+\delta_{jk}R_i
                                   +\delta_{ki}R_j\big)D^2
                                   +R_iR_jR_kD^3\big],
  \end{array}$$ where ${{\sf M}^{}}_{{\rm B}}\equiv{{\sf M}^{0}}_{{\rm B}}$ is the mass of cell B and ${\tilde{\sf M}{}^{2}}_{{\rm B}}\equiv{{\sf M}^{2}}_{{\rm B}}/{{\sf M}^{0}}_{{\rm B}}$ its [*specific*]{} quadrupole moment, while $$\label{Dn}
  D^m\equiv \left.
    \left({1\over r}{\partial\over\partial r}\right)^m
    g(r) \right|_{r=|{\mbox{\boldmath$\scriptstyleR$}}|}.$$ In practice, we calculate the coefficients ${\tilde{{\mbox{\boldmath$\sf C$}}}{}^{m}}_{{{\rm B}}\to{{\rm A}}}
\equiv{{\sf M}^{}}_{{\rm A}}{{\mbox{\boldmath$\sf C$}}^{m}}_{{{\rm B}}\to{{\rm A}}}$, because these obey the symmetry relations ${\tilde{{\mbox{\boldmath$\sf C$}}}{}^{2}}_{{{\rm B}}\to{{\rm A}}}= {\tilde{{\mbox{\boldmath$\sf C$}}}{}^{2}}_{{{\rm A}}\to{{\rm B}}}$ and ${\tilde{{\mbox{\boldmath$\sf C$}}}{}^{3}}_{{{\rm B}}\to{{\rm A}}}=-{\tilde{{\mbox{\boldmath$\sf C$}}}{}^{3}}_{{{\rm A}}\to{{\rm B}}}$ (for $p=3$), which arise from the mutuality of gravity. Since we consider always both directions of any interaction, exploiting these relations substantially reduces the operational effort.

Accumulating Taylor Coefficients {#sec:accum}
--------------------------------

After, for each cell, the coefficients of all its interactions have been accumulated as ${\tilde{{\mbox{\boldmath$\sf C$}}}{}^{m}}_{{\rm A}}= \sum_{{\rm B}}{\tilde{{\mbox{\boldmath$\sf C$}}}{}^{m}}_{{{\rm B}}\to{{\rm A}}}$, where the sum includes all interaction partners B of cell A, we transform back to ${{\mbox{\boldmath$\sf C$}}^{m}}_{{\rm A}}={\tilde{{\mbox{\boldmath$\sf C$}}}{}^{m}}_{{\rm A}}/{{\sf M}^{}}_{{\rm A}}$. Next, for each body, the Taylor series of all relevant cells (those that contain the body), have to be accumulated by first translating to a common expansion center and then adding coefficients.

In contrast to expansions in spherical harmonics, the translation of the Cartesian expansion [(\[local-expn\])]{} to a different center ${{\mbox{\boldmath$z$}}}$ is straightforward. Let ${{\mbox{\boldmath$\sf C$}}^{m,p}}_0$ be the coefficients for expansion center ${{\mbox{\boldmath$z$}}}_0$, then the coefficients for expansion center ${{\mbox{\boldmath$z$}}}_1$ are $$\label{translate}
  {{\mbox{\boldmath$\sf C$}}^{m,p}}_1 = \sum_{n=0}^{p-m} {1\over n!}\,({{\mbox{\boldmath$z$}}}_0-{{\mbox{\boldmath$z$}}}_1)^{(n)} \odot
  {{\mbox{\boldmath$\sf C$}}^{m+n,p}}_0.$$

The Multipole Acceptance Criterion (MAC) {#sec:mac}
----------------------------------------

For the expansion [(\[local-expn\])]{} to converge, we must have $|{{\mbox{\boldmath$x$}}}-{{\mbox{\boldmath$y$}}}-{{\mbox{\boldmath$R$}}}|
<|{{\mbox{\boldmath$R$}}}|$ for all body-body interactions ‘caught’ by a single cell-cell (or cell-body) interaction (see Appendix A). In order to ensure this, we first obtain, for each node, an upper limit $r_{\rm max}$ for the distance of any body within the node from the center of mass (bodies naturally have $r_{\rm max}=0$). We take $r_{\rm max}$ to be either the distance from the cell’s center of mass to its most distant corner [@SW], or [@Benz] $$\max_{{\rm child\;nodes}\;i}\{r_{i,\rm max}+|{{\mbox{\boldmath$z$}}}-{{\mbox{\boldmath$z$}}}_i|\},$$ whatever is smaller. Then $|{{\mbox{\boldmath$x$}}}-{{\mbox{\boldmath$y$}}}-{{\mbox{\boldmath$R$}}}|<\theta|{{\mbox{\boldmath$R$}}}|\;\forall\;{{\mbox{\boldmath$x$}}}\in{{\rm A}}\;{\rm
  and}\;{{\mbox{\boldmath$y$}}}\in{{\rm B}}$, i.e. the nodes are [*well-separated*]{}, if $$\label{well-separated}
  |{{\mbox{\boldmath$R$}}}| > r_{\rm A,crit} + r_{\rm B,crit}
  \quad{\rm with}\quad
  r_{\rm crit} = r_{\rm max} / \theta.$$ The [*tolerance parameter*]{} $\theta$ controls the accuracy of the approximation: for Newtonian forces, the error made in $d$ dimensions by the $p$th order expansion of the form [(\[local-expn\])]{} is ([equation [(\[error:acc\])]{}]{}) $$\begin{aligned}
 \label{error:rel}
  |{\cal R}_p({\mbox{\boldmath$\nabla$}}\Phi_{{{\rm B}}\to{{\rm A}}})| &\le&
  {(p+1)\theta^p\over(1-\theta)^2}\;{{{\sf M}^{}}_{{\rm B}}\over R^2}
\\   \label{error:abs} &\propto&
  {\theta^{p+2}\over(1-\theta)^2}\;r_{\rm B,max}^{d-2}
  \propto
  {\theta^{p+2}\over(1-\theta)^2}\;{{{\sf M}^{}}_{{\rm B}}}^{(d-2)/d}\end{aligned}$$ where ${{\sf M}^{}}_{{\rm B}}\propto r_{\rm B,max}^d$ has been assumed. Since ${{\sf M}^{}}_{{\rm B}}/R^2$ is, to lowest order, the acceleration due to the interaction, [equation [(\[error:rel\])]{}]{} tells us that with constant $\theta$, which is standard tree-code practice, the [*relative*]{} error introduced by every single interaction is approximately constant. [Equation [(\[error:abs\])]{}]{}, however, shows that for the most important case of $d=3$ this practice results in larger [*absolute*]{} errors for interactions with bigger and hence on average more massive cells. This implies that these interactions will dominate the [*total*]{} error of any body’s acceleration. It is, therfore, expedient to balance the absolute errors of the individual interactions, which can be approximately achieved by using a mass-dependent tolerance parameter as follows $$\label{open:mass}
        {\theta^{p+2}\over(1-\theta)^2} =
        {\theta_{\rm min}^{p+2}\over(1-\theta_{\rm min})^2}
        \left({{{\sf M}^{}}\over{{\sf M}^{}}_{\rm tot}}\right)^{(2-d)/d}$$ where ${{\sf M}^{}}_{\rm tot}$ is the mass of the root cell, while $\theta_{\rm min}=
\theta({{\sf M}^{}}_{\rm tot})$ is the new tolerance parameter. If $r_{\rm A,max}\sim
r_{\rm B,max}$ and ${{\sf M}^{}}_{{\rm A}}\sim{{\sf M}^{}}_{{\rm B}}$, this method results in approximately constant absolute acceleration errors. Note that [equation [(\[open:mass\])]{}]{} results in a very weak decrease of $\theta$ with increasing mass: $\theta\propto{{\sf M}^{-1/15}}$ for $d=3$, $p=3$ and $\theta\ll1$.

Instead of [equation [(\[error:abs\])]{}]{}, one can obtain a stricter error limit incorporating the first $p+1$ multipole moments (see Appendix A), which may be turned into an MAC [@WS]. However, our choice [(\[open:mass\])]{} is (i) much simpler and (ii) overcomes already the main disadvantage of $\theta=\,$const, the variations of the absolute individual errors.

Direct Summation {#sec:sum}
----------------

For small $N$ the exact force computation via direct summation is not only more accurate than approximate methods but also more efficient. Therefore, we replace the approximate technique by direct summation whenever the latter results in higher accuracy at the same efficiency, see Appendix B for more details.

If an interaction is executed by direct summation, the Taylor coefficients of the interacting cell(s) are not affected, but the coefficients ${\tilde{{\mbox{\boldmath$\sf C$}}}{}^{0}}$ and ${\tilde{{\mbox{\boldmath$\sf C$}}}{}^{1}}$ of the bodies within the cell(s) are accumulated.

The Algorithm {#sec:alg}
=============

Tree Building and Preparation {#sec:build}
-----------------------------

In the first stage, a hierarchical tree of cubic cells is build, as described in [@BH], albeit cells containing $s$ or less bodies are not divided. Next, the tree is linked such that every cell holds the number of cell children[^2], a pointer to its first cell child, as well as the number of body children, of all body descendants, and a pointer to the first body child. Since cell as well as body children are contiguous in memory (to arrange this we actually use tiny copies of the bodies, called ‘souls’, that only hold a pointer to their body and its mass, position, acceleration and potential) these data allow fast and easy access to all cell and body children of any given cell and to all bodies contained in it.

Next, the masses ${{\sf M}^{}}$, centers of mass ${{\mbox{\boldmath$z$}}}$, radii $r_{\rm max}$ and $r_{\rm
  crit}$, and specific quadrupole moments ${\tilde{\sf M}{}^{2}}$ of each cell are computed in a recursive way from the properties of the child nodes.

A new Generic Tree-Walk Algorithm {#sec:walk}
---------------------------------

One important new feature of our code is the mutual treatment of all interactions: both interacting nodes are source and sink simultaneously. The standard tree-walk, as for instance implemented by the generic algorithm given in [@WS], as well as the the usual FMM coefficient accumulation algorithm, e.g. in [@FMMb], contain an inherent asymmetry between sinks and sources, and thus cannot be used for our purposes. Instead, our algorithm approximates the forces in two steps: an [*interaction phase*]{}, incorporating [equation [(\[coeffs\])]{}]{}, and an [*evaluation phase*]{}, incorporating [equation [(\[local-expn\])]{}]{}.

### The Interaction Phase {#sec:iact}

Because of the mutuality of the interactions, we cannot accumulate the Taylor coefficients ${\tilde{{\mbox{\boldmath$\sf C$}}}{}^{n}}$ ‘on the walk’, but each node must accumulate the coefficients of all its interactions in its own private memory. Cells need storage for ${\tilde{{\mbox{\boldmath$\sf C$}}}{}^{0}}$ to ${\tilde{{\mbox{\boldmath$\sf C$}}}{}^{3}}$, while bodies only need to accumulate ${\tilde{{\mbox{\boldmath$\sf C$}}}{}^{0}}$ and ${\tilde{{\mbox{\boldmath$\sf C$}}}{}^{1}}$, i.e. potential and acceleration. The accumulation of these coefficients is done by the following algorithm with the root cell for arguments $A$ and $B$.

try to perform the interaction between[node]{}s $A$ and $B$;(it cannot be performed)   [**if**]{} ($A=B$)   [**for**]{}(all pairs $\{a,b\}$ of child[node]{}s of $A$)   [Interact]{}($a,b$);   [**else if**]{} ($r_{\rm max}(A) > r_{\rm max}(B)$)   [**for**]{}(all child[node]{}s $a$ of $A$)   [Interact]{}($a,B$);   [**else**]{}   [**for**]{}(all child[node]{}s $b$ of $B$)   [Interact]{}($A,b$);

Thus, if an interaction cannot be executed, using the formulae of the last section – see Appendix B for details – it is split. In case of a mutual interaction, the bigger node is divided and up to eight new mutual interactions are created, while a self interaction of a cell results in up to 36 new interactions between its child nodes. In practice, we use a non-recursive code incorporating a stack of interactions.

### The Evaluation Phase {#sec:eval}

Finally, the Taylor coefficients relevant for each body are accumulated and the expansion is evaluated at every body’s position. After transforming ${\tilde{{\mbox{\boldmath$\sf C$}}}{}^{}}$ to ${{\mbox{\boldmath$\sf C$}}^{}}$ for every cell and body, this is done by the following recursive algorithm, which is initially called with the root cell and an empty Taylor series as arguments.

$T_A$ =[Taylor series]{} due to the ${{\mbox{\boldmath$\sf C$}}^{n}}$ of[cell]{} $A$; translate center of $T_0$ to center of mass of $A$; $T_A$ += $T_0$; (all[body]{} children of $A$) {   evaluate $T_A$ at[body]{}’s position;   add to[body]{}’s potential and acceleration; } [**for**]{}(all[cell]{} children $C$ of $A$)   [Evaluate Gravity]{}($C$, $T_A$);

Thus, the coefficients of the Taylor-series that is eventually evaluated at some body’s position have been added up from all hierarchies of the tree and hence account for all interactions of all cells that contain the body.

Error Assessment {#sec:assess}
================

Two types of force errors are involved in collisionless $N$-body simulations of stellar dynamics. One is the unavoidable error between the smooth force field of the underlying stellar system modeled and the forces [*estimated*]{} from the positions of $N$ bodies (which are sampled from this stellar system). This [ *estimation error*]{} can be reduced by increasing the number $N$ of bodies in conjunction with a careful [*softening*]{}: the Newtonian Greens function $g({{\mbox{\boldmath$x$}}})=G/|{{\mbox{\boldmath$x$}}}|$ is replaced with a non-singular function that approaches the Newtonian form for $|{{\mbox{\boldmath$x$}}}|$ larger than the softening length $\epsilon$. But at fixed $N$, it cannot be decreased below a certain optimum value [@Merr; @Db].

The other type of error is introduced by an approximate rather than exact computation of these estimated forces. While this [*approximation error*]{} can be reduced to (almost) any size (at the price of increasing computational effort), it is sufficient to reduce it well below the level of the estimation error. We will now first assess the approximation error alone and then consider the total error.

Accuracy of the Approximation {#sec:accur}
-----------------------------

We estimate the accuracy of the approximated forces for various choices of the parameters controlling the algorithm and for two typical astrophysical situations: a spherical Hernquist [@Hern] model galaxy, which has density and force per unit mass $$\label{hernquist}
  \rho({{\mbox{\boldmath$x$}}}) = {{\sf M}_{\rm tot} \over 2\pi}
  {r_0\over |{{\mbox{\boldmath$x$}}}| \big(r_0+|{{\mbox{\boldmath$x$}}}|\big)^3},\qquad
  {{\mbox{\boldmath$F$}}}({{\mbox{\boldmath$x$}}})  = - {{{\mbox{\boldmath$x$}}}\over|{{\mbox{\boldmath$x$}}}|}
  { G{\sf M}_{\rm tot}\over\big(r_0+|{{\mbox{\boldmath$x$}}}|\big)^2},$$ and a group of five such galaxies at different positions and with various scale radii $r_0$. In either case, we sample $N=10^4$, $10^5$, and $10^6$ bodies. We use standard Plummer softening, where $g(r)=G/\sqrt{r^2+\epsilon^2}$, with softening lengths $\epsilon$ chosen such as to minimize the estimation error [@Db]. In order to single out the approximation error, we compare with the forces obtained by a computation via direct summation (in double precision; in case of $N=10^6$ for the first $10^5$ bodies only). As measure for the relative error, we compute for each body [@CM] $$\label{error}
  \varepsilon \equiv |a_{\rm approx} - a_{\rm direct}| / a_{\rm direct},$$ where $a$ denotes the magnitude of the acceleration. [Figure \[fig:error\]]{} plots the mean relative error, $\overline\varepsilon$, and that at the 99th percentile, $\varepsilon_{99\%}$, versus the CPU time needed by an ordinary PC (Pentium III/933MHz/Linux/compiler: gcc version 2.95.2) for both constant and mass-dependent tolerance parameter. This figure allows several interesting observations.

1.  For the same $\theta=\,$const and stellar system, the errors decrease with increasing $N$. This is because, at constant [*relative*]{} force error per individual interaction (as is the case for $\theta=\,$const), the total error of some body’s force scales with the inverse square root of the number of individual interactions contributing, which increases with $N$.

2.  At the same operational effort, as measured by the CPU time, the mass-dependent tolerance parameter employing [equation [(\[open:mass\])]{}]{} results in smaller errors than $\theta=\,$const, for $N>10^4$. This advantage becomes more pronounced for larger $N$, because, with $\theta=\,$const, the [ *absolute*]{} force errors of individual interactions contributing to some body’s force vary stronger with increasing $N$ (due to the larger range of cell masses), such that balancing them becomes more beneficial.

3.  Finally, a relative error of $\varepsilon_{99\%}$ of a few per cent or $\overline\varepsilon$ of a few $10^{-3}$ at $N=10^5$, which is commonly accepted to be sufficient in astrophysical contexts [@CM], requires a tolerance parameter $\theta=\,$const$\,\simeq0.65$ or $\theta_{\rm
        min}\simeq0.5$.

The Total Force Error {#sec:error}
---------------------

The important question here is for which choices of $\theta$ is the approximation error negligible compared to the estimation error? To answer this question, we have performed some experiments using samples of $N=10^4$, $10^5$ and $10^6$ bodies drawn from a Hernquist model and computed the mean averaged squared error (MASE) of the force (per unit mass) [@AT]: $$\label{mase}
  {\rm MASE}({{\mbox{\boldmath$F$}}}) = \left\{{1\over{\sf M_{\rm tot}}} \sum_i \mu_i
    \big(\hat{{{\mbox{\boldmath$F$}}}}_i - {{\mbox{\boldmath$F$}}}({{\mbox{\boldmath$x$}}}_i)\big)^2 \right\}.$$ Here, $\hat{{{\mbox{\boldmath$F$}}}}_i$ and ${{\mbox{\boldmath$F$}}}({{\mbox{\boldmath$x$}}})$ are, respectively, the approximately estimated force for the $i$th body and the true force field of the stellar system. The curly brackets denote the ensemble average over many possible random realizations of the same underlying stellar density by $N$ bodies. We used $10^7/N$ ensembles and computed the MASE$({{\mbox{\boldmath$F$}}})$[^3] for various values of $\theta$, but always at optimum $\epsilon$ [@Db]. As one might already have guessed from the behavior of the approximation error in [Figure \[fig:error\]]{}, the relative increase of the MASE$({{\mbox{\boldmath$F$}}})$ is negligible: even for $\theta=0.9$, the approximation error contributes less than one per cent, in agreement with the findings of [@AT]. Based on this result, one may advocate the usage of tolerance parameters larger than $\theta_{\rm
  min}\simeq0.5$. However, the distribution of approximation errors is not normal, and the rms error, which is essentially measured by (the square root of) the MASE$({{\mbox{\boldmath$F$}}})$, may well underestimate the danger of using large $\theta$. We therefore cannot recommend using $\theta_{\rm min}{\mathrel{\hbox{\rlap{\hbox{\lower4pt\hbox{$\sim$}}}\hbox{$>$}}}}0.7$.

Performance Tests {#sec:perform}
=================

Scaling with $N$ {#sec:scaling}
----------------

We measured the CPU time consumption for both constant $\theta=0.65$ and $\theta=\theta({{\sf M}^{}})$ with $\theta_{\rm min}=0.5$, using $s=6$ in either case. [Figure \[fig:N\]]{} plots the consumed time (averaged over many experiments) per body versus $N$ plotted on a logarithmic scale. For the case of $\theta=\theta({{\sf M}^{}})$, [Table \[tab:N\]]{} gives the number of cells as well as the number of individual body-body (B-B), cell-body (C-B) and cell-cell (C-C) interactions, where the latter two are split into those done via a Taylor expansion and direct summation (subscripts ‘app’ and ‘dir’), respectively.

[@rrrrrrrr]{} $N$ & $N_{\rm cells}$ & B-B & C-B$_{\rm app}$ & C-B$_{\rm dir}$ & C-C$_{\rm app}$ & C-C$_{\rm dir}$ & total 1000 & 370 & 210 & 3746 & 2049 & 4057 & 1600 & 11662 3000 & 1090 & 472 & 16870 & 5184 & 30101 & 4472 & 57099 10000 & 3558 & 753 & 61873 & 11914 & 163554 & 11887 & 249981 30000 & 10607 & 494 & 170634 & 24984 & 572322 & 20102 & 790474 100000 & 35282 & 112 & 429035 & 73350 & 1954508 & 53466 & 2516146 300000 & 106065 & 66 & 1041836 & 205821 & 5457445 & 137138 & 6859114 1000000 & 353342 & 1 & 2918326 & 621802 & 16105065 & 393311 & 19023391 3000000 & 1060650 & 2 & 7771584 & 1689115 & 42974890 & 1047627 & 53645379

The tree code requires ${\cal O}(N\log N)$ operations, corresponding to a rising straight line in [Fig. \[fig:N\]]{}. For our code, however, there is a turn-over at $N\sim10^4$, above which the CPU time per body approaches a constant[^4] (for $\theta=\,$const) or even decreases with $N$ (for $\theta=\theta({{\sf M}^{}})$), i.e. the total number of operations becomes ${\cal O}(N)$ or less, which is also evident from the number of interactions in [Table \[tab:N\]]{}.

In order to understand these scalings, let us first consider the simpler case of $\theta=\,$const and a homogeneous distribution of bodies. Then, eight-folding $N$ is equivalent to arranging eight copies of the old root cell into the octants of the new root cell [@BH], and the total number of interactions rises from $N_I$ to $8N_I+N_+$, where $N_+$ interactions are needed for the mutual forces between these octants. In terms of a differential equation, this gives $$\label{ode}
        {{\rm d} N_I\over{\rm d} N}
        \simeq  {N_I\over N} {\Delta\ln N_I\over\Delta\ln N}
        \approx {N_I\over N} + {N_+\over N8\ln 8},$$ where the first term on the right-hand side accounts for the [ *intra-domain*]{} and the second for the [*inter-domain*]{} interactions. [Equation [(\[ode\])]{}]{} has solution $$\label{sol}
        N_I \simeq c_0 N + {N\over8\ln8} \int {N_+\over N^2}\;{\rm d} N.$$ In the tree code, every body requires a constant number of additional interactions, i.e. $N_+\propto N$, and the second term in [(\[sol\])]{} becomes $\propto N\ln N$ dominating $N_I(N)$. However, in the new code, $N_+$ grows sub-linear for large $N$, since a constant number of cell-cell interactions accounts for most new interactions of all bodies. In this case, the second term on the right-hand side of [equation [(\[sol\])]{}]{} also grows sub-linear with $N$ and $N_I$ will eventually be dominated by the first term. That is, in contrast to the tree code, the inter-domain interactions are neglible at large $N$ when compared to the intra-domain interactions. The transition value of $N$ will depend on the tolerance parameter $\theta$ and the distribution of bodies.

Another way of estimating the scaling of the computational costs with $N$ is similar to the FMM approach [@FMMa; @FMMb]: On each level $l$ of the tree there are $\sim8^l$ cells of mass ${{\sf M}^{}}\sim8^{-l}$, i.e. $n({{\sf M}^{}})\propto{\sf
  M}^{-2}$. Each cell has $\propto\theta^{-3}$ interactions, and thus $$\label{scaling}
        N_I \propto \int { d{{\sf M}^{}}\over{\sf M}^{2}\theta^3}.$$ Hence, for $\theta=\,$const, $N_I\propto1/{\sf M}_{\rm min}\propto N$, while a shallower scaling results if $\theta({{\sf M}^{}})$ increases towards smaller masses. Empirically we find for $N>30000$ that the CPU time used by the force computation alone (without tree building) is very well fit by the power-law $\propto N^{0.929\pm0.001}$.

Comparison with Other Methods used in Stellar Dynamics {#sec:comp}
------------------------------------------------------

For various computational techniques used in stellar dynamics, [Figure \[fig:var\]]{} plots $N$ versus the CPU time consumption normalized by $N$, both on logarithmic scales. An ordinary direct method running on a general-purpose computer is slower than our code for any $N{\mathrel{\hbox{\rlap{\hbox{\lower4pt\hbox{$\sim$}}}\hbox{$>$}}}}100$. The [GRAPE-5]{} system [@GRAPE] obtains a $\sim100$ times higher performance by wiring elementary gravity into special-purpose hardware and, for $N{\mathrel{\hbox{\rlap{\hbox{\lower4pt\hbox{$\sim$}}}\hbox{$<$}}}}10^4$, is faster than any other method.

The new code presented here is the fastest method running solely on general-purpose computers and the only one faster than ${\cal O}(N\log N)$. In particular, it out-performs the popular tree code by a factor of 10 and more. At $30000{\mathrel{\hbox{\rlap{\hbox{\lower4pt\hbox{$\sim$}}}\hbox{$<$}}}}N{\mathrel{\hbox{\rlap{\hbox{\lower4pt\hbox{$\sim$}}}\hbox{$<$}}}}3\times10^7$, the only technique that requires less CPU time is a combination of the tree code with a [GRAPE-5]{} board [@GRAPE; @G-tree]. Here, the speed-up due to the usage of special-purpose hardware does not quite reach that for direct methods, because of tree-building and other overheads that cannot be done on the [GRAPE]{} board.

Comparison with Fast Multipole Methods {#sec:fmm}
--------------------------------------

Because our code relies on cell-cell interactions, it may be considered a variant of FMM, introduced by Greengard & Rokhlin [@FMMa; @FMMa1]. However, it differs in several ways from most implementations of FMM. (i) the expansions are centered on the cells’ centers of mass instead of the geometrical centers; (ii) a Cartesian Taylor series is used instead of an expansion in spherical harmonics; (iii) the interaction partners are determined by a multipole acceptance criterion rather than by their mutual grid position; (iv) the mutuality of interactions is fully exploited; (v) the expansion order $p$ is fixed and the accuracy controlled by the parameter $\theta$.

[@rrrlrrl]{} & & & & & & 20000 &13.3 & 233 &$7.9\cdot10^{-4}$ & 0.97 & 136 &$3.7\cdot10^{-4}$ 50000 &27.7 & 1483 &$5.2\cdot10^{-4}$ & 2.64 & 924 &$3.3\cdot10^{-4}$ 200000 &158 & 24330 &$8.4\cdot10^{-4}$ &10.77 & 14694 &$3.4\cdot10^{-4}$ 500000 &268 &138380 &$7.0\cdot10^{-4}$ &29.42 & 91134 &$3.7\cdot10^{-4}$ 1000000 &655 &563900 &$7.1\cdot10^{-4}$ &58.34 &366218 &$3.5\cdot10^{-4}$

To assess the effect of these differences, we compared our code directly with the 3D adaptive FMM code by Cheng, Greengard & Rokhlin [@FMMb]. We performed a test identical to one reported by these authors ($N$ bodies randomly distributed in a cube) on an identical computer (a Sun UltraSPARC with 167MHz) using the same error measure (eq. (57) of [@FMMb]): $$E = \left[\sum_i(\Phi_i - \tilde\Phi_i)^2 \Big/ \sum_i \Phi_i^2\right]^{1/2}$$ where $\Phi$ and $\tilde\Phi$ are the potential computed by direct summation and the approximate method, respectively (unfortunately, Cheng et al. do not give the more relevant error of the accelerations). [Table \[tab:compare\]]{} gives, in the last three columns, the CPU time ($T_{\rm approx}$) in seconds and error $E$ for our code with $\theta=1$ and $s=6$ as well as the time needed for direct summation in 64bit precision ($T_{\rm direct}$); columns 2-4 report the data from Table I of [@FMMb]. On average our code is faster by more than a factor of ten and twice as accurate (even though we compromised the approximation of the potential by omitting the octopole contributions).

What causes this enormous difference? Since for the direct summation the timings are much more similar, we can rule out differences in hardware, compiler, etc.  as cause. An important clue is the fact that our code cannot compete with the Cheng et al. FMM in the regime $E{\mathrel{\hbox{\rlap{\hbox{\lower4pt\hbox{$\sim$}}}\hbox{$<$}}}}10^{-6}$. While our code tries to reach this goal by decreasing $\theta$, FMM obtains it by increasing $p$. In general, the accuracy as well as the performance are controlled by both the order $p$ and the choice of interaction partners, parameterized in our code by $\theta$. Hence, maximal efficiency at given accuracy is obtained at a unique choice of $(p,\theta)$. Apparently, low orders $p$ are optimal for low accuracies, while high accuracies are most efficiently obtained with high orders, instead of an increasing number of interactions (decreasing $\theta$).

Thus, traditional FMM is less useful in the low-accuracy regime, such as needed in stellar dynamics, in agreement with earlier findings [@CM], and our code may be called a variant of FMM optimized for low-accuracy. Clearly, however, a code for which $p$ and $\theta$ can be adapted simultaneously would be superior to both.

Discussion and Summary {#sec:summ}
======================

Our code for the approximate computation of mutual long-range forces between $N$ bodies extends the traditional Barnes & Hut [@BH] tree code by including cell-cell interactions, similar to fast multipole methods (FMM). However, unlike most implementations of FMM, our code is optimized for comparably low-accuracy, which is sufficient in stellar dynamical applications ([§\[sec:fmm\]]{}).

As a unique feature, our code exploits the mutual character of gravity: both nodes of any interaction are sink and source simultaneously. This results in exact conservation of Newton’s third law and substantially reduces the computational effort, but requires a novel tree-walking algorithm ([§\[sec:walk\]]{}) which preserves the natural symmetry of each interaction. Note that the generic algorithm given in [§\[sec:walk\]]{} is not restricted to long-range force approximations, but can be used for any task that incorporates mutuality, for instance neighbor and collision-partner searching.

Complexity {#sec:scal}
----------

The new code requires ${\cal O}(N)$ or less operations ([Fig. \[fig:N\]]{}) for the approximate computation of the forces of $N$ mutually gravitating bodies. For stellar dynamicists, this is the first competitive code better than ${\cal
  O}(N\log N)$. A complexity of ${\cal O}(N)$ was expected for methods based on cell sinks (implying cell-cell interactions) [@FMMa; @FMMb; @WS], but, to my knowledge, hardly ever shown empirically in three dimensions.

Our code obtains a complexity of less than ${\cal O}(N)$ by employing a mass-dependent tolerance parameter $\theta$. The traditional $\theta=\,$constant results in equal [*relative*]{} errors of each interaction, such that the total error of any body’s force is dominated by the interactions with the most massive cells. By slightly increasing the tolerance parameter for less massive cells, we obtain (approximately) equal [*absolute*]{} errors, resulting in a lower total force error than the traditional method for the same number of interactions. Thus, at the same error, we require less interactions. The additional interactions, arising when increasing $N$, occur at ever less massive cells and hence at ever larger tolerance parameters. This causes the computational costs to rise sub-linear with $N{\mathrel{\hbox{\rlap{\hbox{\lower4pt\hbox{$\sim$}}}\hbox{$>$}}}}10^4$ (depending on the accuracy requirements).

Performance {#sec:per}
-----------

We have shown that on general-purpose computers our code out-performs any competitor code commonly used in the field of stellar dynamics. A recent adaptive 3D implementation of FMM [@FMMb] is also out-performed by a factor of ten ([§\[sec:fmm\]]{}), which is related to the fact that traditional FMM codes appear to be good only in the high-accuracy regime. The code presented here was optimized for low accuracies, but by increasing the expansion order $p$, one can easily obtain a version suitable for high accuracies, and it remains to be seen how it would perform compared to traditional FMM.

Currently, the only faster method appears to be a [GRAPE]{}-supported tree code [@G-tree], which uses the special-purpose [GRAPE]{} hardware [@GRAPE]. Unfortunately, unlike the tree code, our code cannot be combined with the current [GRAPE]{} hardware. There are, however, no conceptual obstacles against hard-wiring [equations [(\[hard-coeffs\])]{}]{} into special-purpose hardware, which should yield a speed-up comparable to that of tree to [GRAPE]{} tree, i.e. a factor $\sim50$.

Publication of the Code {#sec:pub}
-----------------------

Our code is written in [C++]{}, also includes a purely two-dimensional version (not described here), and is linkable to [C]{} and [FORTRAN]{} programs. The code will be is publicly available from the author upon request.

A full $N$-body code based on this force algorithm is available under the [ NEMO]{} [@nemo] package (http://bima.astro.umd.edu/nemo). Here, we give an error estimate for the Taylor series approximation of gravity. Using the integral form of the Taylor series remainder, we find for the remainder in [equation [(\[Taylor\])]{}]{}, with ${\mbox{\boldmath$\Delta$}}\equiv({{\mbox{\boldmath$x$}}}-{{\mbox{\boldmath$z$}}}_{{\rm A}})-({{\mbox{\boldmath$y$}}}-{{\mbox{\boldmath$z$}}}_{{\rm B}})=
{{\mbox{\boldmath$x$}}}-{{\mbox{\boldmath$y$}}}-{{\mbox{\boldmath$R$}}}$, $$\begin{aligned}
 \label{remain:pot}
  {\cal R}_p(g) &=&  {{\mbox{\boldmath$\Delta$}}^{(p+1)}\over p!}\,\odot
  \int_0^1 dt\,(1-t)^p\,{\mbox{\boldmath$\nabla$}}^{(p+1)} g({{\mbox{\boldmath$R$}}}+{\mbox{\boldmath$\Delta$}} t)
\\[1ex] \label{remain:acc}
  {\mbox{\boldmath$\nabla$}}{\cal R}_p(g) &=& {{\mbox{\boldmath$\Delta$}}^{(p)}\over(p-1)!}\,\odot
  \int_0^1 dt\,(1-t)^{p-1}\,{\mbox{\boldmath$\nabla$}}^{(p+1)} g({{\mbox{\boldmath$R$}}}+{\mbox{\boldmath$\Delta$}} t),\end{aligned}$$ For Newtonian gravity, $$\begin{aligned}
  \left\| {1\over p!}\int_0^1dt\,(1-t)^p\,
    {\mbox{\boldmath$\nabla$}}^{(p+1)}g({{\mbox{\boldmath$R$}}}+{\mbox{\boldmath$\Delta$}} t) \right\|
  &\le&   {1\over(R-|{\mbox{\boldmath$\Delta$}}|)\,R^{p+1}},          \\[1ex]
  \left\| {1\over (p-1)!}\int_0^1dt\,(1-t)^{p-1}\,
    {\mbox{\boldmath$\nabla$}}^{(p+1)}g({{\mbox{\boldmath$R$}}}+{\mbox{\boldmath$\Delta$}} t) \right\|
  &\le&   {(p+1)R - p|{\mbox{\boldmath$\Delta$}}| \over (R-|{\mbox{\boldmath$\Delta$}}|)^2\,R^{p+1}}.\end{aligned}$$ For the summation over the source cell, we find $$\Big\| \sum_{{\mbox{\boldmath$\scriptstyley$}}_i\in{{\rm B}}} \mu_i {\mbox{\boldmath$\Delta$}}^{(p)} \Big\|
  \le \sum_{k=0}^p {p\choose k} r_{\rm A,max}^k\,\big\|{{\mbox{\boldmath$\sf M$}}^{(p-k)}}_{{\rm B}}\big\|
  \le (r_{\rm A,max}+r_{\rm B,max})^p {{\sf M}^{}}_{{\rm B}}$$ with ${{\sf M}^{}}_{{\rm B}}$ the mass of cell B. With $\theta>(r_{\rm A,max}+r_{\rm
  B,max})/R$, we finally get $$\begin{aligned}
 \label{error:pot}
  |{\cal R}_p(\Phi_{{{\rm B}}\to{{\rm A}}})| &\le&
  {\theta^{p+1}\over1-\theta} {{{\sf M}^{}}_{{\rm B}}\over R} \\ \label{error:acc}
  |{\cal R}_p({\mbox{\boldmath$\nabla$}}\Phi_{{{\rm B}}\to{{\rm A}}})| &\le&
  {(p+1)\theta^p\over(1-\theta)^2}        {{{\sf M}^{}}_{{\rm B}}\over R^2}.\end{aligned}$$

Here, we give the interaction details for [Algorithm 1]{} in [§\[sec:iact\]]{}. Mutual body-body interactions are done by elementary gravity, while body self-interactions are ignored. Mutual interactions between nodes containing $N_1$ and $N_2$ bodies are treated as follows.

1.  If $N_1 N_2 < N_{\rm nn}^{\rm pre}$, execute the interaction by direct summation; otherwise,

2.  if the interaction is well-separated, execute it using [equations [(\[hard-coeffs\])]{}]{}; otherwise,

3.  if $N_1 N_2 < N_{\rm nn}^{\rm post}$, execute the interaction by direct summation; otherwise,

4.  the interaction cannot be executed, but must be split.

Here, $N_{\rm nn}^{\rm pre}$ and $N_{\rm nn}^{\rm post}$ have different values depending whether it is a cell-body (cb) or cell-cell (cc) interaction (see below). Cell self-interactions are done slightly differently:

1.  If $N_1 < N_{\rm cs}$, execute the interaction by direct summation; otherwise,

2.  the interaction cannot be executed, but must be split.

The numbers $N_{\rm cb}^{\rm pre}$, $N_{\rm cb}^{\rm post}$, $N_{\rm cc}^{\rm
  pre}$, $N_{\rm cc}^{\rm post}$, and $N_{\rm cs}$ determine the usage of direct summation. After some experiments, I found the following values to result in most efficient code at given accuracy (but this certainly depends on the implementation details). $$\begin{array}{lcrlcrlcr}
    N_{\rm cb}^{\rm pre}  &=&  3,& 
    N_{\rm cc}^{\rm pre}  &=&  0,& \\
    N_{\rm cb}^{\rm post} &=&128,&
    N_{\rm cc}^{\rm post} &=& 16,&
    N_{\rm cs}            &=& 64.
  \end{array}$$ Note that cell-body interactions with as many as 128 bodies in the cell hardly ever occur, since the interaction algorithm favors interactions between roughly equally sized nodes. Thus, cell-body interactions will almost always be executed.

I thank J. Makino from Tokoy University for helpful discussions and R. Ibata and C. Pichon from Strasbourg Observatory for hardware support.

J.E. Barnes and P. Hut, A hierarchical ${\cal O}(N\log N)$ force calculation algorithm, [*Nature*]{} [**324**]{}, 446 (1986) L. Greengard and V. Rokhlin, A fast algorithm for particle simulations, [*J. Comput. Phys.*]{} [**73**]{}, 325 (1987) L. Greengard and V. Rokhlin, A new version of the fast multipole method for the Laplace equation in three dimensions, [*Acta Numerica*]{} [**6**]{}, 229 (1997) H. Cheng. L. Greengard, and V. Rokhlin, A fast adaptive multipole algorithm in three dimensions, [*J. Comput. Phys.*]{} [ **155**]{}, 468 (1999) R. Capuzzo-Colcetta and P. Miocchi, A comparison between the fast multipole algorithm and the tree code to evaluate gravitational forces in 3D, [*J. Comput. Phys.*]{} [**143**]{}, 29 (1998) W. Dehnen, A very fast and momentum-conserving tree code, [*Astrophys. J.*]{} [**536**]{}, L39 (2000) M.S. Warren, and J.K. Salmon, A portable parallel particle program, [*Comput. Phys. Comm.*]{} [**87**]{}, 266 (1995) J.K. Salmon, and M.S. Warren, Skeletons from the Treecode Closet, [*J. Comput. Phys.*]{} [**111**]{}, 136 (1994) W. Benz, R.L. Bowers, A.Q.W. Cameron, and W.H. Press, Dynamic mass exchange in doubly degenerate binaries I. 0.9 and 1.2 $M_\odot$ stars, [*Astrophys. J.*]{} [**348**]{}, 647 (1990) D. Merritt, Optimal Smoothing for N-body Codes, [*Astron. J.*]{} [**111**]{}, 2462 (1996) W. Dehnen, Towards optimal softening in three-dimensional N-body codes - I. minimizing the force error [*Mon. Not. Royal Astron. Soc.*]{} [**324**]{}, 273 (2001), L. Hernquist, An analytic model for spherical galaxies and bulges, [*Astrophys. J.*]{} [**356**]{}, 359 (1990) E. Athanassoula, E. Fady, J.C. Lambert, and A. Bosma, Optimal softening for force calculations in collisionless N-body simulations, [*Mon. Not. Royal Astron. Soc.*]{} [**314**]{}, 475 (2000) A. Kawai, T. Fukushige, J. Makino, and M. Taiji, GRAPE-5: a special purpose computer for N-body simulations, [*Publ. Astron. Soc. Japan*]{} [**52**]{}, 659 (2000) J. Makino, A tree code with special purpose processor, [ *Publ. Astron. Soc. Japan*]{} [**43**]{}, 621 (1991) P.J. Teuben, The Stellar Dynamics Toolbox NEMO, in Astronomical Data Analysis Software and Systems IV, ed. R. Shaw, H.E. Payne and J.J.E. Hayes, PASP Conf Series [**77**]{}, 398 (1995)

[^1]: One of the most recent members of the FMM family, presented by Cheng, Greengard & Rokhlin [@FMMb], does still not show clear ${\cal O}(N)$ behavior at $N=10^6$ (see, e.g., table I of their paper).

[^2]: Hereafter ‘child’ means a direct sub-node of a cell, while ‘descendant’ refers to any node contained within a cell, including the children, the grand-children and so on.

[^3]: The force field of the Hernquist model has a central singularity, causing a 100% force error at $r=0$, which cannot be resolved by $N$-body methods. In our experiments, we have therefore restricted the summation in [equation [(\[mase\])]{}]{} to $|{{\mbox{\boldmath$x$}}}_i|>\epsilon/2$.

[^4]: The slow rise can be entirely attributed to the tree-building, which is an ${\cal O}(N\log N)$ process.
