The tree and the multipole expansion
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):
\[ \frac{s}{d} < \theta \]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:
Put a cube around all the particles. This is the root.
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.
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
endBuilding 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:
\[ \varphi(\vec d) = -G\sum_k \frac{m_k}{\left|\vec d - \vec s_k\right|} \]Everything now comes from expanding that denominator.
Expanding \(1/|\vec d - \vec s|\)
Write the distance out using the cosine rule:
\[ \left|\vec d - \vec s\right| = \sqrt{d^2 - 2ds\cos\gamma + s^2} \]so
\[ \frac{1}{|\vec d - \vec s|} = \frac{1}{d}\,\frac{1}{\sqrt{1 - 2\left(\frac{s}{d}\right)\cos\gamma + \left(\frac{s}{d}\right)^2}} \]Now comes the moment where the whole thing becomes beautiful. That square root is exactly the generating function of the Legendre polynomials:
\[ \frac{1}{\sqrt{1-2xt+t^2}} = \sum_{n=0}^{\infty}P_n(x)\,t^n, \qquad |t|<1 \]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
\[ A = -2xt + t^2 \]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:
\[ A = -2xt + t^2, \qquad A^2 = 4x^2t^2 + \mathcal{O}(t^3) \]Substituting:
\[ (1+A)^{-1/2} = 1 - \frac{1}{2}\left(-2xt + t^2\right) + \frac{3}{8}\left(4x^2t^2\right) + \mathcal{O}(t^3) = 1 + xt + \left(\frac{3x^2}{2} - \frac{1}{2}\right)t^2 + \mathcal{O}(t^3) \]and reading off the coefficients:
\[ P_0(x) = 1, \qquad P_1(x) = x, \qquad P_2(x) = \frac{3x^2-1}{2} \]which are exactly the first three Legendre polynomials. Setting \(t = s/d\) and \(x = \cos\gamma\):
\[ \boxed{\;\frac{1}{|\vec d - \vec s|} = \frac{1}{d}\sum_{n=0}^{\infty}\left(\frac{s}{d}\right)^n P_n(\cos\gamma)\;} \]Note also the convergence condition \(|t| < 1\), i.e. \(s < d\). The series simply does not converge if the field point is inside the cell. That is not a technicality to wave away. It is why \(\theta\) has to be kept below about 1, and it shows up later as a measured effect.
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:
\[ \varphi(\vec d) = -\frac{G}{d}\left[\underbrace{\sum_k m_k}_{n=0} + \frac{1}{d}\underbrace{\sum_k m_k s_k\cos\gamma_k}_{n=1} + \frac{1}{d^2}\underbrace{\sum_k m_k s_k^2\,\frac{3\cos^2\gamma_k-1}{2}}_{n=2} + \dots\right] \]Let me take these one at a time.
\(n=0\), the monopole. This one is trivial:
\[ \sum_k m_k = M \]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
\[ \sum_k m_k s_k\cos\gamma_k = \hat n\cdot\sum_k m_k\vec s_k \]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\).
\[ \boxed{\text{The dipole term vanishes identically.}} \]\(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\):
\[ s_k^2\,\frac{3\cos^2\gamma_k - 1}{2} = \frac{3(\hat n\cdot\vec s_k)^2 - s_k^2}{2} \]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\):
\[ \frac{3(\hat n_i s_{k,i})(\hat n_j s_{k,j}) - s_k^2}{2} = \frac{1}{2}\hat n_i\hat n_j\left(3s_{k,i}s_{k,j} - s_k^2\delta_{ij}\right) \]Summing over the particles in the cell, everything that depends on the cell's internal arrangement collects into a single object:
\[ \boxed{\;Q_{ij} = \sum_k m_k\left(3\,s_{k,i}\,s_{k,j} - s_k^2\,\delta_{ij}\right)\;} \]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:
| term | what it measures | what you would say |
|---|---|---|
| \(n=0\), monopole | total mass \(M\) | "there is something there, and this is how heavy it is" |
| \(n=1\), dipole | where the mass sits | "and it is over there, not where you said" |
| \(n=2\), quadrupole | is it a cigar or a pancake? | "and it is squashed this way" |
| \(n=3\), octupole | is it lopsided? | "and it is fatter at one end than the other" |
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.
\[ \underbrace{\text{slope of } 1/r}_{\text{cancels: dipole} = 0} \qquad \underbrace{\text{curvature of } 1/r}_{\text{survives: quadrupole} \sim (s/d)^2} \]What the quadrupole tensor means
\(\mathbf{Q}\) is a symmetric \(3\times3\) matrix, so it has six independent components. It is also traceless:
\[ \text{tr}\,\mathbf{Q} = \sum_k m_k\left(3s_k^2 - 3s_k^2\right) = 0 \](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\):
\[ \boxed{\;\varphi(\vec d) = -\frac{GM}{d} - \frac{1}{2}\frac{G}{d^5}\,\vec d\cdot\mathbf{Q}\cdot\vec d\;} \]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:
\[ d = \sqrt{d_ld_l} \;\Longrightarrow\; \partial_k d = \frac{d_k}{d} \;\Longrightarrow\; \partial_k\left(d^{-5}\right) = -5d^{-6}\cdot\frac{d_k}{d} = -\frac{5d_k}{d^7} \]Second, the derivative of the quadratic form. Using \(\partial_k d_i = \delta_{ki}\) and the symmetry \(Q_{ij}=Q_{ji}\):
\[ \partial_k S = \partial_k\left(d_iQ_{ij}d_j\right) = \delta_{ki}Q_{ij}d_j + d_iQ_{ij}\delta_{kj} = Q_{kj}d_j + d_iQ_{ik} = 2Q_{kj}d_j \]The two terms are equal because \(\mathbf{Q}\) is symmetric, which is where the factor of 2 comes from. Now the product rule:
\[ \partial_k\left(\frac{S}{d^5}\right) = \frac{\partial_k S}{d^5} + S\,\partial_k\left(d^{-5}\right) = \frac{2Q_{kj}d_j}{d^5} - \frac{5\,S\,d_k}{d^7} \]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:
\[ a_k^{\text{quad}} = -\partial_k\left(-\frac{G}{2}\frac{S}{d^5}\right) = \frac{G}{2}\left(\frac{2Q_{kj}d_j}{d^5} - \frac{5Sd_k}{d^7}\right) = \frac{G\,Q_{kj}d_j}{d^5} - \frac{5G}{2}\frac{S\,d_k}{d^7} \]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\)),
\[ \boxed{\;\vec a = -\frac{GM}{d^2}\hat r + \frac{G}{d^4}\mathbf{Q}\cdot\hat r - \frac{5G}{2}\left(\hat r\cdot\mathbf{Q}\cdot\hat r\right)\frac{\hat r}{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.
You can also see why \(\theta\) has to be kept below about 1. The expansion parameter is \(s/d\), and the series is only guaranteed to converge for \(s/d < 1\). Push \(\theta\) past that and you are asking a divergent series for an answer, so adding more terms need not help. That is exactly what the paper finds, and what I will show in the results.
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):
\[ \mathbf{Q} = \sum_{l}\mathbf{Q}_l + \sum_{l}m_l\left(3\vec R_l\vec R_l - R_l^2\mathbf{1}\right) \]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
\[ \vec s_k = \vec R_l + \vec y_k \]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:
\[ \sum_{k\in l} m_k\Big[3\left(R_{l,i}+y_{k,i}\right)\left(R_{l,j}+y_{k,j}\right) - \left|\vec R_l + \vec y_k\right|^2\delta_{ij}\Big] \]Expand both pieces. The first:
\[ 3\left(R_{l,i}+y_{k,i}\right)\left(R_{l,j}+y_{k,j}\right) = 3R_{l,i}R_{l,j} + 3R_{l,i}y_{k,j} + 3y_{k,i}R_{l,j} + 3y_{k,i}y_{k,j} \]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\):
\[ -\left(R_l^2 + 2\vec R_l\cdot\vec y_k + y_k^2\right)\delta_{ij} \]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
\[ \underbrace{\mathbf{Q}_l}_{\text{its own internal shape}} \;+\; \underbrace{m_l\left(3R_{l,i}R_{l,j} - R_l^2\delta_{ij}\right)}_{\text{subcell treated as a point mass}} \]Now sum this over the subcells \(l\). The two pieces can be collected separately, which gives
\[ \mathbf{Q} = \sum_{l=1}^{n_{\text{subcell}}}\mathbf{Q}_l + \sum_{l=1}^{n_{\text{subcell}}}m_l\left(3\vec R_l\vec R_l - R_l^2\mathbf{1}\right) \]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)
endWhy 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
\[ s_\ell = \frac{L}{2^\ell} \]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
\[ \frac{s_\ell}{\theta} < d < \frac{2s_\ell}{\theta} \]Now count how many cells of size \(s_\ell\) fit into that shell. The shell's volume is
\[ V_\ell = \frac{4\pi}{3}\left[\left(\frac{2s_\ell}{\theta}\right)^3 - \left(\frac{s_\ell}{\theta}\right)^3\right] = \frac{4\pi}{3}\cdot\frac{7s_\ell^3}{\theta^3} \]and each cell occupies volume \(s_\ell^3\), so the number of them is
\[ n_\ell \approx \frac{V_\ell}{s_\ell^3} = \frac{28\pi}{3\theta^3} \]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
\[ n_{\text{terms}} \sim \frac{28\pi}{3\theta^3}\log_8 N = C(\theta)\log N \quad\Longrightarrow\quad \text{total cost} \sim N\log N \]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\):
| \(N\) | \(\langle n_\text{terms}\rangle\) | change |
|---|---|---|
| 1024 | 124.0 | |
| 2048 | 150.8 | +26.8 |
| 4096 | 171.5 | +20.7 |
| 8192 | 191.5 | +20.0 |
| 16384 | 205.3 | +13.8 |
| 32768 | 217.6 | +12.3 |
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.
What one particle actually sees
The average hides a lot of structure. Below is the full distribution of \(n_\text{terms}\) across all 32768 particles, for the Plummer model and for a uniform sphere with the same \(N\). This is the paper's Fig. 5:
Left to right is how many terms one particle had to add up, and the height is how many of the 32768 particles needed that many. The peak is narrow, so almost everybody does a similar amount of work, and the tail on the left is the halo particles who get off cheap.
Two things stand out, and both are physics rather than numerics.
The uniform sphere gives a tight, symmetric distribution centred on 122. Every particle sits in a similar neighbourhood, so every particle does a similar amount of work.
The Plummer model is centred higher (218) and has a long tail towards smaller values. The reason becomes obvious when you plot \(n_\text{terms}\) against radius:
Distance from the centre goes right, work done goes up, one dot per particle. Particles in the crowded middle do the most work and particles in the empty halo do the least, which is exactly where the effort should go.
The reason is the geometry overhead. A halo particle has a shallow tree above it and the whole cluster sitting far away, so a handful of big cells covers everything. A core particle has many levels of subdivision right next to it, and every one of those levels has to be opened.
My numbers here are \(217.6\) and \(122.2\); the paper reports \(221\) and \(121\). Given that this is a completely independent implementation from a written description, I am happy with that.
One trap: a particle pulling on itself
There is a subtle bug the paper warns about in section II. It is invisible unless you go looking for it.
Consider a particle sitting near the edge of a large cell that it is itself inside. The cell's centre of mass could be at the far edge, so \(d\) is large-ish and \(s/d < \theta\) might be satisfied. For \(\theta \gtrsim 1\) this really can happen. The walk then accepts the cell as a single lump, but the particle's own mass is part of that lump. The particle exerts a force on itself.
The fix is a geometric check: if the particle lies inside the cell, always open it, regardless of \(s/d\).
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])
endI have this on by default. The paper investigates it and finds the effect negligible for \(\theta \lesssim 1.2\), which matches what I see, but it costs almost nothing to be correct.
The walk, in full
Putting it all together, here is the complete force calculation for one particle, the piece of code that dominates the entire runtime:
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
endNote also the softening in the quadrupole terms. The paper applies an ad hoc substitution \(r^4 \to (r^2+\varepsilon^2)^2\); I instead use \(r^n \to (r^2+\varepsilon^2)^{n/2}\) throughout, which has the advantage of being exactly \(-\nabla\) of a consistently softened quadrupole potential. The two agree when \(\varepsilon = 0\) and differ only in how they behave at large softening, which is one of the things I test.
Where we are
We now have the complete method:
Build an octree over the particles, \(O(N\log N)\), done once per step.
Fill in mass, centre of mass and quadrupole for every cell by a bottom-up sweep, \(O(N)\).
For each particle, walk the tree accepting cells with \(s/d < \theta\), \(O(\log N)\) each.
Push everything forward with leapfrog.
The next part is the payoff: does it actually behave the way the counting argument says, and how wrong are the forces?
Previous: Part 2, moving the particles
Next: Part 4, results
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.