Worlds from Dust — Algorythmic Explorations

Part one built a structure for the far field: a pyramid of summed-up mass, because in a galaxy distant stars only matter in bulk. Part three built one for the near field: a sorted grid of buckets, because a drop of water only feels what's touching it. Nothing stops one simulation from using both, with gravity from the pyramid and contacts from the grid. That gets you dust that attracts, collides and sticks, which is roughly how planets start.

The ring at the top of the page is 20,000 dust grains on nearly circular orbits around a star. Every grain pulls on every other grain through part one's pyramid, and pushes on the grains it's actually touching through part three's sort. Left alone, the ring grows clumps, and the clumps collide and grow into moonlets. It's a simplified 2D version of how planetesimals form.

Getting there means seeing why gravity on its own can't build anything solid, getting part three's grid out of its box, and running both structures over the same particles every step.

Step 1 · No surfacesGravity alone can't build a planet

Part one introduced softening as a fix for dividing by zero. It also has a physical meaning: it says bodies have no surface. Inside the softening length, gravity fades out instead of blowing up, and two particles aimed straight at each other pass through like ghosts. For a galaxy that's correct. Stars are so small compared to the distances between them that when two galaxies collide, almost no stars actually hit each other. Astronomers call systems like that collisionless.

Planets don't form that way. Here's a cold, slowly spinning disc of 20,000 dust grains with gravity and nothing else:

The cloud falls inward, overshoots, and swells back out. Gravity is a conservative force, so whatever energy the fall gains, the rebound hands back. A cloud with nothing but gravity has no way to get rid of the energy of its own collapse, so however long you leave it, it stays a loose cloud and never becomes anything solid. Press both to turn on contacts and the same cloud does something a galaxy can't. That needs part three's grid, and first the grid has to get out of its box.

Step 2 · Open spaceA grid with no walls

Part three's grid depended on a fixed box. Each cell's index came from its position inside the box, and the box is what kept the number of cells finite. The water tank had walls anyway, so that cost nothing. A collapsing dust cloud is in open space, though, and an orbiting ring throws the odd grain a long way out. Covering that with a dense grid means choosing between a huge box, with millions of counters for empty space, and a tight box that particles escape from.

A spatial hash keeps the cells and drops the box. A cell is still a pair of integers, which column and which row, but now anywhere in the infinite plane. Those two integers get hashed into a fixed-size table of buckets:

fn bucketOf(c: vec2i) -> u32 {
  let h = (u32(c.x) * 0x9E3779B1u) ^ (u32(c.y) * 0x85EBCA77u);
  return h & (TABLE - 1u);   // TABLE is a power of two
}

At the default 20,000 grains the table has 65,536 buckets, and the demo switches to 262,144 and then 1,048,576 as you add grains. Nothing else in part three's pipeline changes. The histogram counts particles per bucket, the scan turns the counts into start positions, and the scatter groups particles by bucket.

The cost is that unrelated cells, possibly far apart, sometimes land in the same bucket, so the 3×3 neighbour search picks up impostors. That doesn't matter, because the force pass already measures the distance to every candidate. A real neighbour passes the distance test and a grain from the far side of the ring doesn't, so the check that was already there throws impostors out. The one thing to watch is that two of the nine cells can hash to the same bucket. The kernel removes duplicate buckets before searching, or it would count those contacts twice.

The demo below makes the impostors visible. It's part three's pile of 3,500 discs, but its 24×24 = 576 cells share a hash table that starts at a deliberately tiny 128 buckets. Yellow is the 3×3 window around your cursor. Red cells are unrelated cells that happen to share one of its buckets, and red discs are the impostors the distance test throws away. Slide the table size up and the red gets rare.

Step 3 · One bufferTwo structures, one buffer

Putting the two together was mostly plumbing. Part three's sort reorders the particles every substep, so grain #4,071 is a different grain after each shuffle. Part one's solver kept each body's mass in a separate array, indexed by body number. Bolt the two together without thinking and the masses stay put while the particles move, so after every substep each grain has some other grain's mass.

The usual fix is to sort the masses along with the particles (a key–value sort, one more array to shuffle). I went with something simpler: every grain weighs the same. Dust doesn't need the variety a galaxy does. If all the masses are equal, reordering the particles can't change the gravity, so the pyramid just counts grains per cell and the mass buffer goes away. The total mass of the dust is fixed, and the grain-count slider only changes how finely it's divided up.

With that, both structures run every substep over the same buffer:

sort     → count, scan ×3, scatter                     // part three's hash sort: 5 dispatches
pyramid  → clear, bounds, scatter, resolve, reduce ×8  // part one's pyramid: 12 dispatches
force    → tree walk + 9 cells, integrate              // new: 1 dispatch

That's 18 dispatches a substep and four substeps a frame. Here's what the force kernel sees, one structure at a time:

The two structures aren't the same size, and they shouldn't be. The pyramid is 256×256 at its finest level and refits its root box to the cloud every substep, because gravity needs to see everything. The hash cells stay one grain diameter wide, because contacts happen at the scale of a grain. The force kernel is the only place they meet: walk the pyramid for the pull of the whole cloud, check the nine cells around the grain for anything it's touching, add the two, and integrate. The force pass runs on the freshly sorted buffer, so a grain's contact neighbours are also its neighbours in memory.

The star isn't in the buffer at all. It's one analytic GM/r² term in the force kernel, the same way part two's dome was a force term rather than something stored with the particles. There's also a soft leash past a fixed radius that pulls stray grains back. Without it, one grain thrown far out would stretch the pyramid's root box and make every cell coarser.

Step 4 · DissipationCollisions that lose energy

The data structures don't supply the last piece. Here's the collapsing cloud again. Try all three buttons:

Gravity only is the ghost cloud from Step 1. Contacts only is a handful of marbles in space: the cloud spins apart, because nothing holds it together. Both collapses and stays collapsed. Grains pile onto grains and the swirl tightens into a dense, slowly spinning clump. That's a planetesimal.

The difference comes from one term in the contact force. The contact model is the same as part three's sand, a spring–dashpot. The spring pushes overlapping grains apart. The dashpot is a damper, like a piston pushed through oil, and it resists the two grains moving toward or away from each other:

acc += nrm * (dia - r) * P.stiffness;   // spring: push apart by the overlap
let vn = dot(vel - parts[k].zw, nrm);   // how fast the pair closes or separates
acc -= nrm * vn * P.damping;            // dashpot: take some of that away

The dashpot makes every collision inelastic. Each bump turns a little of the ordered falling motion into nothing (in real dust it would become heat). The stickiness slider is the damping coefficient, from 0 to 80 with 30 as the default. Since the dashpot also resists grains moving apart, turning it up makes them tend to stay stuck together. Slide it to zero and the collisions become elastic: gravity plus perfectly bouncy contacts still conserves energy, so the collapse turns into a bouncing swarm instead of a clump. Stars in a galaxy have no way to shed energy like this, which is why galaxies stay puffy. Dust does, and that's what lets it settle into something solid.

Aside · Newton's pendulums and Safronov's book

Newton measured this property for the Principia. He swung pendulum balls of different materials into each other and recorded how much of their closing speed survived the bounce: about 5/9 for tightly wound wool, about 15/16 for glass. That ratio is now called the coefficient of restitution, and turning up stickiness turns it down.

The theory that planets grow this way, with dust settling, colliding and sticking and the bigger bodies sweeping up the rest, was worked out largely by Viktor Safronov in Moscow in the 1950s and 60s. He collected it in a Russian monograph in 1969. A NASA translation followed in 1972, and that's how most Western planetary scientists came to know it. It became the basis of the standard picture of planetesimal growth, and simulations like this one still follow its loop.

Step 5 · The diskA ring of dust

The collapsing cloud is planet formation on easy mode, since the only way to go is in. A more realistic nursery is a rotating disk, where angular momentum resists collapse. Here's the ring from the top of the page with the controls showing. The grains slider goes from 10,000 to 150,000.

What fights clumping in a disk is shear. Inner grains orbit faster than outer ones, so any clump that starts to form gets smeared into an arc as its inner edge pulls ahead of its outer edge. Watch the first few seconds: the ring doesn't clump evenly. It grows tight spiral streaks where self-gravity pulls material together while shear stretches it out. Where enough grains crowd together, a knot condenses out of a streak, collisions damp its internal motion, and it becomes a moonlet. From then on it sweeps its lane, scattering or absorbing what it meets.

Aside · Saturn's rings do this

Planetary scientists call the late stage of this oligarchic growth, a name from Eiichiro Kokubo and Shigeru Ida's 1998 simulations. Those followed up to 4,000 planetesimals, on special-purpose hardware from the same Tokyo project as part one's GRAPE machines. The slider above goes to 150,000.

Saturn's A ring runs a version of the experiment where we can watch. It's a sheet of colliding ice particles whose self-gravity keeps trying to build clumps while shear keeps tearing them apart. Cassini saw both outcomes: self-gravity wakes, short-lived elongated clumps like the streaks above, and moonlets about a hundred metres across that stay intact and stir up propeller-shaped disturbances in the ring around them.

This leaves a lot out. Real planet formation is 3D. Real dust sits in gas, which drains orbits through drag and, through the streaming instability, herds grains into clumps. Real collisions can shatter grains as well as stick them. The dust here is also much heavier relative to its star than in any real disk (in the code's units the dust totals 0.12 and the star 0.066), which makes self-gravity far stronger than it would be. The skeleton matches the research codes, though: a tree for the far field, a sorted grid for the near field, and collisions that lose energy. Each missing piece is one more force or one more sorted array.

EpilogueWhat carries over

Structures compose. The pyramid and the hash grid never talk to each other. Each answers one question, "what does everything far away add up to?" or "what am I touching?", and the force kernel asks both.
Dissipation is what lets things settle. Gravity conserves energy, so a cloud with nothing else can clump but can't become solid. One damping term in the contact force is enough to turn clouds into planetesimals and rings into moonlets.
Change the data so the structures can't disagree. The sort reorders particles, and gravity used to look masses up by index. Making every grain identical removed the problem instead of managing it, the same kind of move as part one's pointer-free tree.

Things I'd try next: a third coordinate, which both structures take the same way part two's octree did, to turn the ring into a 3D disk; letting colliding grains merge, which needs particles to disappear and masses to vary, so stream compaction and a key–value radix sort (both built on scans); gas drag, which is one more force; and longer timesteps or multirate integration, since real accretion takes an enormous number of orbits and this page squeezes it into seconds.

Next: Part five, Gravity in Frequency Space. Gravity solved with an FFT, Poisson's equation reduced to one division per wave, and a cosmic web grown from small ripples.