3. The tree and the multipole expansion

Part 3 of my recreation of Hernquist (1987). This is the long one, and the one that actually contains the physics of the method. Everything so far was setup.

The question: how do I compute the force on \(N\) particles without doing \(N^2\) work?

The idea, in one picture

Suppose I want the force on one particular star, and there is a globular cluster of ten thousand stars sitting a long way away. Do I need all ten thousand terms?

Obviously not. From far enough away the cluster is just a lump of mass \(M\) at a particular place. One term instead of ten thousand.

The question is: how far is far enough? The answer has to be relative, because a cluster of size \(s\) seen from distance \(d\) looks compact if \(s \ll d\). So the natural thing to compare is the ratio \(s/d\), which is exactly the paper's eq. (1.1):

Notice that \(s/d\) is roughly the angle the cell subtends on the sky. So \(\theta\) is literally "how big is a clump allowed to look before I stop treating it as a point".

So, a picture. A cell of width \(s\) sitting at distance \(d\) covers an angle of roughly \(s/d\) on your sky. The test asks whether that angle is small enough to ignore what is inside it.

One phrase in that picture needs unpacking, since it refers to something built in the next section. Every cube of space in this method gets cut into 8 equal smaller cubes, in the same way a Rubik's cube is eight little cubes stacked two by two by two. Those eight are called the cube's children. "Look at the 8 children" therefore means: this clump is too close and too wide to be trusted as a single blob, so throw the blob away, take the eight smaller cubes it is made of, and ask the same question of each one separately. Some of those will be far enough to accept; the ones that are not get chopped into eight again. In two dimensions the pictures below use 4 children instead of 8, for the obvious reason that a square splits into four.

Notice the criterion is an angle, not a distance. A clump twice as wide may be treated as one lump only if it is also twice as far away. That is exactly how vision works, and it is why the right thing to test is \(s/d\) rather than \(s\) or \(d\) on its own.

To turn this into an algorithm I need two things: a way to organise the particles into nested clumps of every size (the tree), and a way to describe a clump by more than just its total mass (the multipole expansion). Let me do the tree first, because it is the easy half.

Building the tree

The Barnes-Hut recipe is a recursive subdivision of space:

1. Put a cube around all the particles. This is the root.

2. If a cell has more than one particle in it, cut it into 8 equal sub-cubes (in 3D, hence octree) and hand each particle to whichever sub-cube contains it.

3. Repeat until every cell has at most one particle.

The result is a hierarchy of cells: one enormous cell containing everything, eight cells of half the size, sixty-four of a quarter and so on down to individual particles. Cells in dense regions get subdivided many times; cells in empty regions stop early.

In 2D that looks like this (a quadtree, four children per cell instead of eight, same logic and easier to look at). 500 particles drawn from the Plummer sphere:

Dots are the particles, boxes are the cells the tree made. Notice nobody told it where the crowd is: it worked that out on its own.

The subdivision is deep and fine in the dense middle, shallow out in the sparse halo. This is what makes tree codes so much more flexible than grid methods, which have to commit to a resolution in advance.

That picture is the finished tree though, and it hides the thing that confused me at first, which is that the tree is not designed, it is grown. So here is the same construction with the particles going in one at a time:

Red is the particle being inserted right now, blue are the ones already placed. Watch what happens when a red dot lands in a box that is already occupied: the box immediately splits into four, and both particles fall into their own quarters. Sometimes they land in the same quarter again and it has to split again, which is why one insertion can suddenly add several levels. The title tracks the count, and by the end 120 particles have produced 212 cells and 8 levels. Nobody chose that depth, it is just where the splitting stopped.

The Julia is a straightforward insert-one-particle-at-a-time loop. The one subtlety is the splitting you just watched: when a particle arrives at a cell that already holds one, the sitting particle has to be pushed down a level first:

function insert!(t::Octree, pos::AbstractMatrix, i::Integer)
    node = 1
    while true
        if t.leafpart[node] == 0 && all(t.child[k, node] == 0 for k in 1:8)
            t.leafpart[node] = i          # empty cell: park it here
            return
        elseif t.leafpart[node] != 0
            # occupied leaf: push the sitting tenant down one level first
            j = t.leafpart[node]
            t.leafpart[node] = 0
            oj = octant(t, node, pos[1, j], pos[2, j], pos[3, j])
            c = newnode!(t, child_center(t, node, oj)..., 0.5 * t.size[node])
            t.child[oj, node] = c
            t.leafpart[c] = j
            # ... and carry on with particle i in the same cell
        else
            o = octant(t, node, pos[1,i], pos[2,i], pos[3,i])
            if t.child[o, node] == 0
                c = newnode!(t, child_center(t, node, o)..., 0.5 * t.size[node])
                t.child[o, node] = c
                t.leafpart[c] = i
                return
            end
            node = t.child[o, node]
        end
    end
end

Building the tree costs \(O(N\log N)\), and in practice it is a small fraction of the total time. Hernquist quotes under 10%, and I find the same.

Walking the tree

With the tree built, computing the force on particle \(i\) is a walk from the root:

• If the cell is empty, skip it.

• If the cell holds a single particle, add the ordinary pairwise force (and skip it if it is particle \(i\) itself).

• Otherwise compute \(s/d\). If \(s/d < \theta\), accept the cell: add one term for the whole thing and do not look inside. If not, open it and put its eight children on the stack.

That walk, made visible. This is the same 500-particle quadtree, and the star marks the particle we are computing the force on. Orange boxes are cells that got swallowed whole; red dots are particles that had to be summed one by one. Drag the slider to change \(\theta\):

Move the slider to change \(\theta\) and watch the count above the plot. Small \(\theta\) means many small boxes and a lot of work; large \(\theta\) means a few big boxes and very little work.

At \(\theta = 0.2\) the walk is fussy: 161 terms, descending nearly to individual particles even quite far away. At \(\theta = 1.5\) it is down to 15 terms, with one enormous box covering the whole far side of the cluster.

Also notice the spatial pattern. Nearby particles are always handled individually, and the cells used get bigger the further away they are. The method automatically spends its effort where the force is largest and the geometry matters most.

The multipole expansion

Now the real content. When we "accept" a cell, what exactly do we replace it with?

The crude answer is: a point mass at its centre of mass. That is the monopole approximation, and it is what the original Barnes-Hut paper used. We can do better, and doing better is most of what Hernquist's paper is about. This is the part sir wanted us to actually understand, so nothing below is skipped.

The setup

Take one cell. Put the origin at its centre of mass. The cell contains particles of mass \(m_k\) at positions \(\vec s_k\), which are small, of order the cell size \(s\). We want the potential at a field point \(\vec d\) which is far away (\(d \gg s\)). Let \(\gamma_k\) be the angle between \(\vec d\) and \(\vec s_k\).

The exact potential is just the sum over particles:

Everything now comes from expanding that denominator.

Expanding \(1/|\vec d - \vec s|\)

Write the distance out using the cosine rule:

so

Now comes the moment where the whole thing becomes beautiful. That square root is exactly the generating function of the Legendre polynomials:

The first few terms come out by hand in four lines, and doing that made the structure concrete for me. Write the thing under the root as \(1 + A\) with

and use the binomial series \((1+A)^{-1/2} = 1 - \tfrac{1}{2}A + \tfrac{3}{8}A^2 - \tfrac{5}{16}A^3 + \dots\). We want everything up to \(t^2\), so we need \(A\) in full and \(A^2\) only to leading order:

Substituting:

and reading off the coefficients:

which are exactly the first three Legendre polynomials. Setting \(t = s/d\) and \(x = \cos\gamma\):

Note also the convergence condition \(|t| < 1\), i.e. \(s < d\). The series simply does not converge if the field point is too close. That is not a technicality to wave away, and it is worth being careful about which distance has to be small, because the answer is not the one I assumed for a long time. The condition is that \(|\vec s_\alpha|/d < 1\) for every particle \(\alpha\) in the cell, not for some typical size of the cell. I come back to this below, in the section on what the opening test actually tests, and it turns out to change the safe range of \(\theta\) by a factor of nearly two.

I checked the general term in Mathematica rather than trusting my memory beyond \(n=2\):

$Assumptions = d > 0 && s > 0 && -1 <= ct <= 1;
Table[FullSimplify[SeriesCoefficient[1/Sqrt[d^2 - 2 d s ct + s^2], {s, 0, n}]
                   - LegendreP[n, ct]/d^(n + 1)], {n, 0, 4}]
(* -> {0, 0, 0, 0, 0} *)

Term by term

Substituting back into the potential:

Let me take these one at a time.

\(n=0\), the monopole. This one is trivial:

the total mass of the cell. This gives \(-GM/d\), i.e. the cell treated as a point. Nothing about the shape of the cell survives here.

\(n=1\), the dipole. Write \(\hat n = \vec d/d\) for the direction to the field point. Then \(s_k\cos\gamma_k = \hat n\cdot\vec s_k\), so

And now the payoff for a decision made way back at the start. We put the origin at the centre of mass, which is defined by precisely \(\sum_k m_k\vec s_k = 0\).

\(n=2\), the quadrupole. The first term that actually survives, so I am going slowly here. Again using \(s_k\cos\gamma_k = \hat n\cdot\vec s_k\):

Write it in index notation, with \(\hat n\cdot\vec s_k = \hat n_i s_{k,i}\) (summing over repeated indices) and using \(\hat n_i\hat n_j\delta_{ij} = |\hat n|^2 = 1\):

Summing over the particles in the cell, everything that depends on the cell's internal arrangement collects into a single object:

which is exactly the paper's eq. (2.3). The \(n=2\) term is then \(\tfrac{1}{2}\hat n\cdot\mathbf{Q}\cdot\hat n\).

I verified this collapse symbolically too, since it is easy to drop a factor:

nh = {n1, n2, n3}; xk = {a1, a2, a3}; s2 = xk . xk;
legTerm = s2 LegendreP[2, (xk . nh)/Sqrt[s2]];
Qk = Table[mk (3 xk[[i]] xk[[j]] - s2 KroneckerDelta[i, j]), {i, 3}, {j, 3}];
FullSimplify[(mk legTerm - (1/2) nh . Qk . nh) /. n3 -> Sqrt[1 - n1^2 - n2^2]]
(* -> 0 *)

What these terms actually mean

Before pushing on, stop and ask what we have just built. The algebra hides a very simple picture.

The multipole expansion is a way of describing a lump of matter by progressively finer features, in the same way you might describe a person from further and further away:

Each term is a finer detail than the last, and each matters less the further away you stand. That is precisely what the factors of \((s/d)^n\) are saying: detail costs distance. From far enough away, everything is a point.

You have met this ladder before, probably several times without it being labelled:

• The Earth's gravity field. The Earth bulges at the equator, and that bulge is exactly a quadrupole. Satellite people call its coefficient \(J_2\), and every GPS orbit accounts for it.

• Tides. The Moon pulls harder on the near side of the Earth than the far side. Subtract the average and what is left is the quadrupole field.

• The CMB. Decomposing the microwave sky into \(\ell = 0, 1, 2, \dots\) is this same Legendre expansion on a sphere. \(\ell=1\) is the dipole from our own motion.

• Nuclear physics. Deformed nuclei are catalogued by their electric quadrupole moment: same tensor, different force.

Why the error is second order, intuitively

The \((s/d)^2\) law is easier to believe as a picture than as a derivation.

Treating a cell as a point mass at its centre gets the average distance right but nothing else. Now ask what you got wrong. Particles on the near side of the cell are closer than you assumed, so they pull more than your point-mass estimate; particles on the far side are further, so they pull less.

To first order in the cell size those two errors are equal and opposite, so they cancel, because the centre of mass sits exactly at the balance point. That is the dipole vanishing, seen from the other side.

What does not cancel is that \(1/r\) is curved. The extra pull you gain by moving a bit closer is bigger than the pull you lose by moving the same bit further away. So the near side wins by a little, and that residue, the curvature of \(1/r\) rather than its slope, is the quadrupole.

What the quadrupole tensor means

\(\mathbf{Q}\) is a symmetric \(3\times3\) matrix, so it has six independent components. It is also traceless:

(the \(\delta_{ii} = 3\) in three dimensions), which knocks it down to five.

Physically, \(\mathbf{Q}\) measures how far from spherical the cell is. The cleanest way to see it: if the particles in a cell are distributed with perfect spherical symmetry, then by symmetry \(\sum_k m_k s_{k,i}s_{k,j} = \tfrac{1}{3}\delta_{ij}\sum_k m_k s_k^2\), and substituting gives \(Q_{ij} = 0\) exactly.

The result

Putting the surviving terms together, and using \(\hat n\cdot\mathbf{Q}\cdot\hat n = \vec d\cdot\mathbf{Q}\cdot\vec d / d^2\):

which is the paper's eq. (2.2). The acceleration follows by \(\vec a = -\nabla\varphi\). I did that gradient in index notation, since it is exactly the sort of thing I get wrong by a factor of two.

We need \(\partial_k\) of \(S/d^5\), where \(S \equiv d_iQ_{ij}d_j\). Two ingredients. First, the derivative of the distance itself:

Second, the derivative of the quadratic form. Using \(\partial_k d_i = \delta_{ki}\) and the symmetry \(Q_{ij}=Q_{ji}\):

The two terms are equal because \(\mathbf{Q}\) is symmetric, which is where the factor of 2 comes from. Now the product rule:

Therefore, remembering the \(-\tfrac12 G\) prefactor in \(\varphi\) and the overall minus in \(\vec a = -\nabla\varphi\), the two signs cancel and the \(\tfrac12\) eats the 2:

The monopole part is the familiar \(-GMd_k/d^3\). Putting them together and pulling out unit vectors with \(\hat r = \vec d/d\) (so \(Q_{kj}d_j/d^5 = Q_{kj}\hat r_j/d^4\) and \(Sd_k/d^7 = (\hat r\cdot\mathbf{Q}\cdot\hat r)\hat r_k/d^4\)),

which is the paper's eq. (2.4). I did not want to take on trust that this is minus the gradient of the potential above, so:

Qm = {{q11, q12, q13}, {q12, q22, q23}, {q13, q23, q33}};
rv = {x, y, z}; rr = Sqrt[x^2 + y^2 + z^2];
phi = -Mc/rr - (1/2) (rv . Qm . rv)/rr^5;
accelPaper = -Mc rv/rr^3 + (Qm . rv)/rr^5 - (5/2) (rv . Qm . rv) rv/rr^7;
FullSimplify[-Grad[phi, {x, y, z}] - accelPaper]
(* -> {0, 0, 0} *)

Good. The paper's equations are consistent.

How good is it?

The expansion is a power series in \(s/d\). Since the dipole vanishes, the leading term we throw away is:

• monopole only: we drop the \(n=2\) term, so the error is \(\mathcal{O}\!\left((s/d)^2\right)\)

• through quadrupole: we drop the \(n=3\) term, so the error is \(\mathcal{O}\!\left((s/d)^3\right)\)

That is a prediction, so let me test it. I took a cloud of 200 particles packed into a cube of side \(s\), computed the acceleration it produces at distance \(d\) exactly, and compared against both approximations, averaging over 400 directions:

Going left means the clump is further away, going down means the answer is more correct. Blue treats the clump as a single point, orange also adds its shape. Orange sits below blue everywhere and the gap keeps widening to the left, so adding the shape always helps and helps more the further you are. The slopes come out 2 and 3, which is what we derived above, so the maths and the computer agree.

So the plot is really a measurement of two numbers, the two slopes. The orange line is steeper than the blue one, meaning its error dies away faster as you back off, and the gap between the lines is how much the quadrupole buys you at that distance. At the far left, where \(s/d = 0.02\), that gap is a factor of about a hundred.

Fitting slopes on the small-\(s/d\) end gives 2.00 for the monopole and 3.03 for the quadrupole. Those are the \(2\) and \(3\) derived above, and nothing was fitted to make them come out. I was pleased with this one. It is a clean case of theory predicting a number and the computer producing it.

The test does not test what you think it tests

I used to finish the section above with a sentence like this: the expansion parameter is \(s/d\), the series converges only for \(s/d<1\), so keep \(\theta\) below about 1. It sounds fine. It is wrong, and the way it is wrong is worth a section of its own, because it is the single most instructive mistake I made in this whole project.

Go back to the convergence condition. The Legendre series converges when \(|t|<1\) with \(t = |\vec s_\alpha|/d\), and that has to hold for every particle \(\alpha\) in the cell, not for some average or typical size. So the quantity that decides convergence is

the distance from the expansion centre, which is the cell's centre of mass, out to its own furthest particle. The series converges when \(b_{\max}/d < 1\).

But the opening test never looks at \(b_{\max}\). It looks at \(s\), the width of the cube, which is a completely different number. Nothing in \(s/d<\theta\) says anything directly about where the particles inside the cell actually are.

How far apart can they be?

Take the friendly case first. If the centre of mass sits at the geometric centre of the cube, the furthest point of the cube is a corner, at half the body diagonal:

Now the unfriendly case. Nothing pins the centre of mass to the middle. Put nearly all the mass in one corner, so the centre of mass sits essentially at that corner, and leave one straggler particle in the opposite corner. Now \(b_{\max}\) is the full body diagonal:

So the test \(s/d<\theta\) buys us only

and if we want the guarantee \(b_{\max}/d<1\) to follow from the test, we need

Here is the whole argument in one picture. Cut a plane through the cube along a face diagonal, so that the body diagonal lies inside the plane and is drawn at its true length. The section is a rectangle with sides \(s\) and \(\sqrt2\,s\), so its own diagonal is \(\sqrt{1+2}\,s = \sqrt3\,s\): the worst possible \(b_{\max}\). The dashed circle of that radius is where the series stops converging. The three stars sit at \(d = s/\theta\), which is the closest the test will let a field point come, for three values of \(\theta\):

The middle star lands exactly on the circle, and that is the entire content of \(\theta \le 1/\sqrt3\). For anything larger, the acceptance test is willing to put the field point inside the radius of convergence.

So why does it work anyway?

Because the worst case is a caricature, and I can measure how far real cells are from it. I took the \(N = 32768\) tree, walked it with 200 target particles, and recorded \(b_{\max}\) for every cell that was actually accepted.

The typical accepted cell has \(b_{\max}/s = 0.67\) at \(\theta=0.5\) and \(0.74\) at \(\theta=1\). The whole distribution piles up just below \(\sqrt3/2 = 0.866\), which is exactly what you would get from a cell with its centre of mass near the middle and a particle out near a corner: the friendly case. In the entire measurement, across every cell accepted by 200 different particles, not one cell reached \(b_{\max}/s = 1.25\), never mind the worst case of \(1.73\).

The expansion parameter itself tells the same story more directly. Plotting \(b_{\max}/d\) over the accepted cells:

• the worst-case bound \(\sqrt3\,\theta\) reaches 1 at \(\theta = 0.577\), as it must, by construction;

• the largest \(b_{\max}/d\) the tree actually produces anywhere climbs far more slowly, and only reaches 1 near \(\theta \approx 0.83\);

• the ordinary cell is nowhere near trouble at all. Its \(b_{\max}/d\) is still only \(0.38\) at \(\theta = 1\).

And \(0.83\) earns its keep. In part 6 I build a deliberately awkward configuration, a big cluster with a small compact satellite sliding across cell boundaries, and ask when particles start getting badly wrong forces. Nothing at all goes wrong up to \(\theta = 0.7\). The first bad particles appear near \(\theta = 0.9\), just above the measured crossing, and by \(\theta = 1.3\) the worst particle is wrong by \(242\%\).

What each term is actually seeing

Before leaving the expansion, it is worth asking what the terms mean operationally, because that tells you when each one is needed and when it is wasted. The cleanest way I found is an experiment you can do in twenty lines.

Build three cells that the monopole cannot possibly tell apart: same total mass, same centre of mass, same cell. Then look at how wrong each approximation is.

• (a) round. A roughly spherical cloud. By Newton's shell theorem this should pull almost exactly like a point mass, so there is nearly nothing for the higher terms to do.

• (b) a bar. Stretched along one axis. This is a pure quadrupole shape: it is symmetric under \(\vec s \to -\vec s\), so its third moment vanishes.

• (c) lopsided. A heavy clump on one side and a light one on the other, placed so that the centre of mass is still dead centre. This is the shape a quadrupole cannot describe.

The measured force error at \(s/d = 0.5\), averaged over 240 directions:

Read the bottom two rows together, because that is where the lesson is.

The bar and the lopsided cell have almost the same quadrupole, \(0.124\) against \(0.121\), and almost the same monopole error, \(2.4\%\) against \(2.3\%\). Yet the quadrupole fixes the bar by a factor of 28 and the lopsided cell by a factor of only 8. Adding the octupole then does nothing more for the bar (a factor \(1.2\)) and another factor of 3 for the lopsided one.

This is the same ladder as in electrostatics, and it has the same structure everywhere it appears: each term describes a finer feature of the shape than the one before, and each term is blind to whatever first appears at the next order. The monopole cannot see position, the dipole cannot see elongation, the quadrupole cannot see lopsidedness. Since every new term costs memory and arithmetic, where to stop is a question of economics and not of principle, and part 6 answers it with a cost-accuracy frontier.

Building \(\mathbf{Q}\) for every cell, cheaply

There is a practical problem left. Every cell in the tree needs its own \(\mathbf{Q}\), and computing each one directly from its particles would cost \(O(N\log N)\) per level. Wasteful.

The fix is a parallel-axis theorem for the quadrupole. A parent cell's \(\mathbf{Q}\) can be assembled from its children's, and the paper gives it as eq. (2.5):

where \(l\) runs over the subcells, \(m_l\) is the subcell's mass and \(\vec R_l\) is the offset of subcell \(l\)'s centre of mass from the parent's.

Deriving it is worth the five minutes, because the cancellation that makes it work is the same one that killed the dipole. Take a particle \(k\) in subcell \(l\). Its position relative to the parent's centre of mass is

where \(\vec y_k\) is measured from the subcell's own centre of mass. Substituting into the definition of \(\mathbf{Q}\) and summing over the particles of that one subcell:

Expand both pieces. The first:

and the second, using \(|\vec R_l+\vec y_k|^2 = R_l^2 + 2\vec R_l\cdot\vec y_k + y_k^2\):

Now sort the seven terms by how many powers of \(\vec y_k\) they carry, and sum over \(k\in l\):

• No \(\vec y\): \(\;\sum_k m_k\left(3R_{l,i}R_{l,j} - R_l^2\delta_{ij}\right) = m_l\left(3R_{l,i}R_{l,j} - R_l^2\delta_{ij}\right)\), since \(\sum_{k\in l}m_k = m_l\).

• Linear in \(\vec y\): \(\;3R_{l,i}\sum_k m_ky_{k,j} + 3R_{l,j}\sum_k m_ky_{k,i} - 2R_{l,m}\delta_{ij}\sum_k m_ky_{k,m}\). Every one of these contains \(\sum_{k\in l}m_k\vec y_k\), which is zero by the definition of the subcell's centre of mass. All three vanish.

• Quadratic in \(\vec y\): \(\;\sum_k m_k\left(3y_{k,i}y_{k,j} - y_k^2\delta_{ij}\right)\), which is precisely \(\mathbf{Q}_l\), the subcell's own quadrupole about its own centre of mass.

So the whole subcell contributes

Now sum this over the subcells \(l\). The two pieces can be collected separately, which gives

and that is the paper's eq. (2.5), character for character.

It is also exactly the parallel-axis theorem you met for moments of inertia in mechanics, \(I = I_{\text{cm}} + Md^2\): a "shape about its own centre" piece plus a "point mass displaced" piece. Same algebra, same reason.

The shift is not optional, and it took me a moment to see why. Every cell needs its moments about its own centre of mass, because that is what made the dipole vanish. But a parent's centre of mass is not any of its children's, so a child's stored numbers are simply the wrong numbers for the parent and have to be translated first. The alternative is to recompute each cell's moments from its own particles, which costs \(O(n)\) for a cell of \(n\) particles and \(O(N\log N)\) over the whole tree, as expensive as the force calculation it was meant to accelerate. With the shifting rule a parent is built from its eight children in constant time, so the whole tree's moments are filled in with a single bottom-up sweep costing \(O(N)\) in total:

for o in 1:8
    c = t.child[o, node]
    c == 0 && continue
    mc = t.mass[c]
    Rx = t.com[1, c] - t.com[1, node]
    Ry = t.com[2, c] - t.com[2, node]
    Rz = t.com[3, c] - t.com[3, node]
    R2 = Rx*Rx + Ry*Ry + Rz*Rz
    # the child's own shape ...
    qxx += t.quad[1, c]; qxy += t.quad[2, c]; qxz += t.quad[3, c]
    qyy += t.quad[4, c]; qyz += t.quad[5, c]; qzz += t.quad[6, c]
    # ... plus the shape it acquires by sitting off-centre
    qxx += mc * (3Rx*Rx - R2);  qxy += mc * (3Rx*Ry)
    qxz += mc * (3Rx*Rz);       qyy += mc * (3Ry*Ry - R2)
    qyz += mc * (3Ry*Rz);       qzz += mc * (3Rz*Rz - R2)
end

Why this gives \(N\log N\)

Stand on a rooftop and look at a city. The houses on your own street you see one by one; a few streets away you see blocks; further out, neighbourhoods; on the horizon, one smudge. Now count how many things you are looking at in each of those bands. It is roughly the same number in each, because things further away are bigger but subtend the same angle.

Now the counting argument. It came out less obvious than I expected, so I am doing it slowly.

Fix one particle and walk down the tree level by level. At level \(\ell\), cells have size

where \(L\) is the size of the root cell. A cell at that level is accepted when \(s_\ell/d < \theta\), i.e. when its distance satisfies \(d > s_\ell/\theta\), and it is opened otherwise. So the cells actually used at level \(\ell\) are those which are accepted while their parent was not. Parents live at level \(\ell-1\) with size \(2s_\ell\), and were opened when \(d < 2s_\ell/\theta\). The cells contributing terms at level \(\ell\) therefore sit in the spherical shell

Now count how many cells of size \(s_\ell\) fit into that shell. The shell's volume is

and each cell occupies volume \(s_\ell^3\), so the number of them is

The \(s_\ell\) has cancelled completely. Every level of the tree contributes the same number of terms, controlled only by \(\theta\). Which is the crux: doubling the resolution adds one more level, and one more level costs a constant.

The tree has \(\sim\log_8 N\) levels, so

That is the entire promise of the method, and it makes a sharp prediction I can test: doubling \(N\) should add a constant amount to \(n_{\text{terms}}\), not double it. My measurement, at \(\theta = 1\):

Compare with the direct sum, where the same column would read 1023, 2047, 4095, 8191, 16383, 32767. Thirty-two times more particles cost me 1.75 times more work per particle. That is the whole game.

A better count

So the shell argument gets the growth with \(N\) right and the dependence on \(\theta\) badly wrong. That bothered me for a long time and I want to fix it properly, because a factor of two in the cost of a simulation is not a rounding error, and because "real cells are not uniform, so what do you expect" is not an explanation, it is an excuse.

Go back and look at what the shell count quietly assumed. The cluster is infinite, and every cell is full. Neither is true, and each failure bites at a different end of the range.

• A cell has to be inside the cluster. The accepted cells of width \(s\) live in a shell at distances \(s/\theta\) to \(2s/\theta\). For small \(\theta\), and for the big cells near the top of the tree, that shell is pushed out beyond the edge of the cluster, where there is nothing to count. The volume is there; the matter is not. This kills the largest cells.

• A cell has to contain particles. Down at the bottom of the tree, cells are much smaller than the mean spacing between particles, so most of them are empty, and an empty cell contributes nothing and costs nothing. This kills the smallest cells.

The crude count includes both populations at full strength. That is why it overestimates, and why it overestimates worst at small \(\theta\), where the shells are furthest out.

Putting both cut-offs in

Write \(W(d)\) for the average volume of the system lying within a distance \(d\) of one of its own particles. This is the quantity that knows about the edge: for small \(d\) it is just \(\tfrac43\pi d^3\), and once \(d\) is bigger than the system it saturates at the system's whole volume. Then the shell around a target particle contains

cells rather than \(28\pi/3\theta^3\) of them, and of those, the fraction that actually hold anything is

with \(\rho\) the number density: the Poisson probability that a box of volume \(s^3\) is not empty. Summing over the levels, with \(s_\ell = L/2^\ell\),

For a uniform ball of radius \(R\) the overlap volume is known in closed form, and I got it out of a computer algebra system rather than trusting an integral I did by hand:

opened = d2 <= 0 || s * s >= theta2 * d2
if !opened && forced_subdivision
    opened = inside_cell(t, node, pos[1, i], pos[2, i], pos[3, i])
end

function tree_accel(t::Octree, pos, i, theta, eps; quadrupole = false,
                    forced_subdivision = true, stack = Int32[])
    empty!(stack); push!(stack, Int32(1))
    ax = ay = az = 0.0
    nterms = 0
    eps2 = eps * eps; theta2 = theta * theta

    while !isempty(stack)
        node = pop!(stack)
        m = t.mass[node]
        m == 0 && continue

        dx = pos[1,i] - t.com[1,node]
        dy = pos[2,i] - t.com[2,node]
        dz = pos[3,i] - t.com[3,node]
        d2 = dx*dx + dy*dy + dz*dz

        p = t.leafpart[node]
        if p != 0
            p == i && continue           # a particle feels no force from itself
            r2 = d2 + eps2
            f = m / (r2 * sqrt(r2))
            ax -= f*dx; ay -= f*dy; az -= f*dz
            nterms += 1
            continue
        end

        s = t.size[node]
        opened = d2 <= 0 || s*s >= theta2 * d2
        if !opened && forced_subdivision
            opened = inside_cell(t, node, pos[1,i], pos[2,i], pos[3,i])
        end
        if opened
            for o in 1:8
                c = t.child[o, node]
                c != 0 && push!(stack, c)
            end
            continue
        end

        # far enough away: replace the cell by its multipoles
        r2 = d2 + eps2; r = sqrt(r2)
        f = m / (r2 * r)
        ax -= f*dx; ay -= f*dy; az -= f*dz

        if quadrupole
            qxx = t.quad[1,node]; qxy = t.quad[2,node]; qxz = t.quad[3,node]
            qyy = t.quad[4,node]; qyz = t.quad[5,node]; qzz = t.quad[6,node]
            qdx = qxx*dx + qxy*dy + qxz*dz      # Q . d
            qdy = qxy*dx + qyy*dy + qyz*dz
            qdz = qxz*dx + qyz*dy + qzz*dz
            dQd = dx*qdx + dy*qdy + dz*qdz      # d . Q . d
            r5 = r2*r2*r
            c1 = 1 / r5
            c2 = 2.5 * dQd / (r5 * r2)
            ax += c1*qdx - c2*dx
            ay += c1*qdy - c2*dy
            az += c1*qdz - c2*dz
        end
        nterms += 1
    end
    return ax, ay, az, nterms
end