Results
Results
Part 4 of my recreation of Hernquist (1987). The code is built and checked; now let me measure it and see whether the paper's conclusions hold up.
Everything below uses the paper's own parameters: a Plummer model with \(M=1\), \(r_0=0.2\), cutoff \(R=1\), \(G=1\), softening equal to the mean interparticle separation at the half-mass radius and \(\delta t = 0.025\).
Does it stay put?
The first and most basic test. The model was constructed to be in equilibrium, so if the code is right, the cluster should sit there and do essentially nothing for a long time.
\(N = 4096\), \(\theta = 1\), 1000 leapfrog steps. Press play:
Press play and watch whether the blob stays the same size. It does, so the cluster is in equilibrium and the code is not quietly blowing it up.
It breathes a little and jiggles, but it does not collapse and it does not evaporate. That is the paper's Fig. 1, and a harder test than it looks. A bug in the force calculation, the tree or the sampling would show up here immediately as a cluster that explodes or implodes within a few crossing times.
A sharper version of the same test is to track the radii containing fixed fractions of the mass. If the model really is stationary, these should be flat lines:
Time goes right, radius goes up. Each curve is the radius holding a fixed share of the mass, so flat lines mean the cluster is not expanding or collapsing. They are flat.
The small wobble at early times is the one we predicted in part 1 from the cutoff at \(R=1\). The 90% radius is the noisiest, which makes sense: it is set by a handful of particles way out in the sparse halo, so it is the most affected by small-number statistics.
How fast is it?
Now the question the paper is named after.
Cost against \(N\)
Time for one force evaluation, per particle, against \(\log_{10}N\), for a range of \(\theta\). Straight lines on these axes mean \(\text{cost} \propto N\log N\):
For \(\theta \gtrsim 0.4\) the curves are straight, so the \(N\log N\) scaling is real. For \(\theta \lesssim 0.3\) (dashed) they bend upwards, because so few cells are being accepted that the method is drifting back towards \(O(N^2)\). In the limit \(\theta \to 0\) it is a direct sum, by construction.
This reproduces the paper's Fig. 3, including where the behaviour changes. Hernquist concludes that \(\theta \gtrsim 0.4\) is needed to stay in the \(N\log N\) regime, and I find the same.
Cost against \(\theta\)
The same data sliced the other way, which is the paper's Fig. 4:
\(\theta\) goes right, time per force evaluation goes up. The cost drops very fast up to about \(\theta = 1\) and then flattens, so there is no point pushing \(\theta\) higher: you pay in error and get almost no speed.
That flattening is the practical reason \(\theta \sim 0.5\) to \(1\) is the range everybody actually uses.
The direct cause is visible if you plot the number of terms instead of the time:
Same story with the machine taken out of it. \(\theta\) goes right, terms per particle goes up. Thousands of terms at small \(\theta\), about a hundred at large \(\theta\).
Past \(\theta \approx 1\) almost everything is already clumped, so relaxing further just regroups cells that were being accepted anyway.
When is it actually worth it?
A tree code is not free. It has to build a tree, chase pointers and do a comparison at every node, while a direct sum is a tight, perfectly predictable inner loop that a compiler can vectorise beautifully. For small \(N\), the simple method wins.
Number of particles goes right, time goes up. Below the crossing point the simple \(N^2\) sum is actually faster, so a tree is only worth building once you have a few thousand particles.
The crossover at \(\theta = 1\) lands at \(N \approx 1700\), which is exactly what the paper estimates ("the two are comparable for \(N \approx 1700\)"). At \(N = 32768\) my tree is \(24\times\) faster than the direct sum, and the gap keeps widening.
What does the quadrupole cost?
How much extra time the shape term costs, as a ratio. It is about 1.3 to 1.5 times, and it does not grow with \(N\), so the shape term is cheap.
Keeping the quadrupole moments makes the force evaluation roughly 1.3 to 1.5 times more expensive, agreeing with the paper's "a factor \(\sim 1.5\)". Whether that is worth paying is the next question.
How wrong are the forces?
This is the part I care about most. Speed is worthless if the answer is rubbish.
The measure is the paper's, from eqs. (3.4) and (3.5). For each Cartesian component, compare the tree acceleration against a direct sum on the same particles, take the mean offset, then take the mean absolute deviation about it:
\[ \overline{\delta a_i} = \frac{1}{N}\sum_j\left(a^{\text{tree}}_{i,j}-a^{\text{direct}}_{i,j}\right), \qquad A(\delta a_i) = \frac{1}{N}\sum_j\left|a^{\text{tree}}_{i,j}-a^{\text{direct}}_{i,j}-\overline{\delta a_i}\right| \]and quote \(A(\delta a_i)/\bar a_i\) as a percentage.
\(\theta\) goes right, force error goes up. This is the plot that decides everything: about \(0.2\%\) at \(\theta = 0.5\), \(1.5\%\) at \(\theta = 1\), and useless beyond \(\theta = 2\). The usable range is narrow.
This is the paper's Fig. 6, and if you remember one plot from the series, make it this one.
My numbers next to the paper's:
| \(N\) | \(\theta = 0.5\) | \(\theta = 1.0\) |
|---|---|---|
| 1024 | 0.43% | 2.91% |
| 4096 | 0.31% | 2.32% |
| 16384 | 0.22% | 1.71% |
| 32768 | 0.18% | 1.55% |
| paper, \(N=32768\) | \(\approx 0.2\%\) | \(\approx 1.4\%\) |
The error gets better with more particles
Look at the columns above again. The error goes down as \(N\) goes up. That surprised me at first. More particles means a deeper tree and more approximations, so surely more error?
Particles go right, error goes up, both log. The line slopes gently downwards, so adding particles does reduce the error, but so slowly that it is not a practical way to get accuracy.
Over this range the fitted behaviour is \(\sim N^{-1/5}\), exactly the scaling the paper reports (my exponent is \(-0.18\)).
The reasoning I eventually settled on goes like this. Each accepted cell contributes a small error of random sign, so the total error on a particle is a random walk over the \(n_\text{terms}\) cells it used. A random walk of \(n\) steps of size \(\epsilon\) accumulates \(\sqrt{n}\,\epsilon\), not \(n\epsilon\). Meanwhile the true acceleration itself grows as more mass is resolved. The competition between those two gives a slow decline, and it is genuinely slow: to halve the error you would need to increase \(N\) by a factor of 32.
The practical conclusion is unchanged, only sharpened: you cannot fix a tree code by throwing particles at it. You fix it with \(\theta\).
Does the quadrupole help?
Two curves, with and without the shape term. The gap between them is what the shape term buys, and it is large on the left and small on the right. So the shape term is worth most when you are already being careful.
At small \(\theta\), enormously. At \(\theta = 0.5\) the quadrupole version has error \(0.05\%\) against the monopole's \(0.18\%\), nearly four times better. At \(\theta = 1\), \(1.04\%\) against \(1.55\%\), only about a third better.
The trend is clearer as a ratio, which is the paper's Fig. 7:
The same gap drawn as a ratio. It shrinks towards zero as \(\theta\) grows, meaning at large \(\theta\) the extra term stops helping. That is the series failing to converge, exactly as part 3 warned.
This is exactly the convergence issue from part 3. The expansion is a power series in \(s/d\); once \(\theta \gtrsim 1\) we are evaluating it at \(s/d \gtrsim 1\), where it is not guaranteed to converge, so adding the next term is not guaranteed to help.
Softening fights the expansion
The last error test, and the one with the nicest physics in it. What happens when the softening length \(\varepsilon\) becomes comparable to the interparticle separation \(\lambda\)?
Softening goes right, error goes up. A little softening helps, too much hurts, and past \(\varepsilon \approx \lambda\) the shape term actually makes things worse than not using it at all.
For small \(\varepsilon\) the error drops, because softening smooths out the very close pairs whose forces are hardest to approximate. But past \(\varepsilon/\lambda \approx 1\) it climbs again, and strikingly the quadrupole version becomes worse than the monopole one.
This is not a numerical accident. The multipole expansion I derived in part 3 is an expansion for point masses: it assumes each particle in the cell is a point at \(\vec s_k\). Softening replaces every particle by a smeared blob of size \(\varepsilon\). Once \(\varepsilon\) is comparable to the cell size, the "shape" that \(\mathbf{Q}\) is carefully describing is no longer the shape of the actual mass distribution, so the quadrupole is correcting for a geometry that the force law does not have. So the correction stops being a correction.
This is the paper's Fig. 8, and its conclusion that the useful range is \(\varepsilon/\lambda \approx 1\) to \(2\) is what I find too.
The conservation laws
Approximating forces has consequences beyond a number being slightly wrong, and this is the part of the paper I find most instructive.
Energy
Time goes right, energy error goes up. The important thing is that even the exact \(N^2\) sum does not sit at zero, so some of the error was never the tree's fault.
Energy drift turns out to be a genuinely noisy quantity, since the same run with a different random seed can differ by a factor of three, so single numbers are not worth much here. Averaging over eight independent realisations:
| direct sum | \(\theta = 0.5\) | \(\theta = 1\) | |
|---|---|---|---|
| mine, 8 seeds | \(0.22 \pm 0.10\%\) | \(0.26 \pm 0.10\%\) | \(0.96 \pm 0.14\%\) |
| Hernquist (1987) | \(0.20\%\) | \(0.32\%\) | \(0.68\%\) |
The first two land on the paper's values. The \(\theta=1\) column is about two standard deviations high, in the same direction as my slightly larger force error at that \(\theta\). See the discussion on the front page.
The instructive part is that the direct sum does not conserve energy either. Its \(0.12\%\) comes entirely from the leapfrog's finite step size, and it is the floor that no force approximation can go below. At \(\theta = 0.5\) the tree's contribution is comparable to the integrator's; at \(\theta = 1\) the tree dominates. That tells you how to spend effort: there is no point using a fancier integrator at \(\theta=1\).
But that floor is not fixed. It depends on the softening, and much more violently than I expected. The same 800-step run, repeated across a factor of 64 in \(\varepsilon\):
Softening goes right, energy error goes up, both on log axes. Read it right to left: as the softening shrinks the error climbs by a factor of a thousand. So the left end of this plot is not more accurate physics, it is a broken simulation.
| \(\varepsilon/\lambda\) | 0.0625 | 0.125 | 0.25 | 0.5 | 1.0 | 2.0 | 4.0 |
|---|---|---|---|---|---|---|---|
| \(\Delta E/E\) | \(+51.6\%\) | \(+34.0\%\) | \(+15.9\%\) | \(+3.35\%\) | \(+0.15\%\) | \(+0.005\%\) | \(-0.29\%\) |
At \(\varepsilon = \lambda/16\) the simulation gains half its own energy. The mechanism is explained in part 1: \(\varepsilon\) sets the hardest encounter in the system, so a fixed \(\delta t\) can only resolve encounters that are not too hard.
Momentum, the one that really breaks
Now the subtle failure, and the paper's sharpest observation.
Think about what happens when a particle interacts with a distant cell. The particle feels the pull of the whole cell. But the cell does not feel the pull of the particle. Instead, each particle inside that cell separately does its own tree walk, and its view of our particle may be as part of a different cell entirely.
So the force of \(A\) on \(B\) is not equal and opposite to the force of \(B\) on \(A\). Newton's third law is broken by construction. And momentum conservation is a direct consequence of the third law, so momentum is not conserved:
Time goes right, distance wandered goes up. The centre of mass should never move at all, and it does. Lowering \(\theta\) is the only thing that helps.
It is a random walk, and the size of it depends only on \(\theta\). At \(\theta = 1\) the centre wanders about \(1\%\) of the system radius over 1000 steps; at \(\theta = 0.5\) it is roughly ten times smaller.
This is the sort of thing that makes me like this paper. It is not just "here is a fast algorithm", it is "here is precisely which piece of physics I traded away and here is how much".
Is the tree code more collisional?
The last question, and it is a genuinely worrying one.
A star cluster relaxes because individual stars deflect each other. The granularity of the mass distribution slowly randomises orbits, on the relaxation time \(t_r\). A tree code replaces groups of stars with single massive pseudo-particles. Bigger lumps mean bigger deflections. So does the tree make the simulation artificially more collisional than it should be, silently corrupting the physics?
The paper tests this with the Standish & Aksnes method, and so did I: fire massless test particles through the system, measure how much they are deflected in one crossing, then form
\[ t_r = \frac{\langle \Delta t\rangle}{\langle\sin^2\Phi\rangle} \]Then quote everything relative to a direct calculation, \(t_r(\theta)/t_r(0)\):
\(\theta\) goes right, and up is the relaxation time compared with the exact calculation. A value of 1 means the tree changed nothing. It stays at 1 up to about \(\theta = 1.1\), so grouping particles does not secretly make the physics grainier.
The answer is no, at least in the range anyone uses. Only beyond \(\theta \approx 1.2\) does the ratio start to fall, reaching about \(0.8\) at \(\theta = 1.4\) with softening.
This matches the paper's Fig. 10, including the detail that the \(\varepsilon = 0\) curves are much noisier. Without softening, the occasional very close encounter dominates the deflection statistics and the variance goes through the roof.
What I take away from it
Pulling the threads together, the picture that emerges is clean:
The \(N\log N\) scaling is real for \(\theta \gtrsim 0.4\), and the tree beats a direct sum from about \(N \sim 1700\) upwards.
The force error is controllable and predictable, about \(1\%\) at \(\theta=1\) and \(0.2\%\) at \(\theta=0.5\), improving only glacially with \(N\).
The quadrupole is best understood as a way of affording a larger \(\theta\), not as free accuracy.
Energy conservation is comparable to a direct calculation. Momentum conservation is not, and cannot be fixed except by lowering \(\theta\).
Collisionality is unaffected for \(\theta \lesssim 1.1\).
Which lands on \(\theta \approx 0.5\) to \(1.0\) as the sensible operating range, exactly where the paper's abstract puts it.
Some things I did not do, if you want to take this further:
Individual time steps. Every particle currently takes the same \(\delta t\), set by the fastest one. Giving core particles smaller steps than halo particles should be a large win.
Octupole terms. The \(n=3\) term gives \(\mathcal{O}((s/d)^4)\). Worth it, or is the extra bookkeeping self-defeating? I went and did this one; see part 6.
A better opening criterion. \(s/d<\theta\) is crude, since it ignores where the mass inside a cell actually sits. Modern codes use criteria based on the multipole moments themselves.
Reversing the drift. Since momentum is not conserved, why not re-centre the system every step? The paper suggests it; does it introduce artefacts of its own?
Previous: Part 3, the tree and the multipole expansion
Next: Part 5, the code, then Part 6, doing better than 1987
Hope this helps you in some way. If you like it then share with others if possible.
If you have some queries, do let me know in the comments or contact me using the informations that are given on the page About Me.