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 