Gravity in Frequency Space

The simulation at the top of this page is a different kind of thing from the rest of the series. The dots aren't stars, or even galaxies. Each one is a sample of the matter in a box that stands in for a patch of the universe, 589,824 samples in all. The box starts out almost smooth. Gravity grows the small ripples into knots, filaments form between the knots, and the filaments wall off empty voids. Astronomers call this the cosmic web. It's the largest pattern in how matter is laid out, and galaxy surveys see the same structure in the real sky.

None of the earlier methods suit this problem. Part one's tree works by summarizing clumps, and the early universe doesn't have any yet. Part three's grid assumes forces with a cutoff distance, and gravity has none. So this part uses a third structure, a mesh plus a Fourier transform, and solves gravity in frequency space.

Step 1 · Periodic boundariesA universe with no edges

A galaxy is a thing in empty space. You can draw a box around it and the box ends. The universe doesn't end: wherever you stand, it carries on in every direction, statistically the same. Simulations handle that by making the box periodic. Leave through the right wall and you come back in through the left, and a particle near one edge pulls on particles near the opposite edge through the wall. Part four freed the grid from its box. Cosmology glues the box back together into a torus.

That's awkward for part one's tree. On a torus every particle also feels every periodic copy of the box, forever, and that sum converges badly. Tree codes handle it with Ewald summation, a 1921 technique from crystallography, at a real cost in time and code. There's a second problem. In the early universe the density is nearly uniform, so the force on each particle is the small amount left over after big pulls from every direction cancel. The tree's approximation errors are small next to the big pulls, but not next to that leftover.

What works better is to stop thinking about gravity as pairs and treat it as a field. Mass density ρ is a function over the box. The gravitational potential φ is another function, and Poisson's equation relates them:

∇²φ = 4πG ρ        // the potential's curvature is set by the density

∇² is the Laplacian, the sum of the second derivatives in each direction. In 2D it's the field's east–west curvature plus its north–south curvature, which measures how much the value at a point sags below or bulges above its neighbours.

The force is then the downhill direction, F = −∇φ. There are no pairs. Each particle reads the slope of a field that already includes every mass in the box and every periodic copy of it. The job is to get the potential everywhere from the density everywhere. On a 512×512 mesh that's 262,144 equations, each unknown tied to its neighbours, and the cosmology demo solves them twice per frame. That sounds like a lot of work. In frequency space the core of the solve is three lines.

Step 2 · FourierAny shape is a sum of waves

Any function over the box, however jagged, can be written as a sum of sine waves. That's Fourier's theorem, and the demo below lets you try it. It takes a signal of 256 samples and rebuilds it from its waves, lowest frequency first. The bump only needs a handful. The square wave needs far more, and however many you add, the rebuilt version overshoots right next to the jumps. Noise comes back too once the slider is near the top: 256 samples and 128 waves (each with an amplitude and a phase) are 256 numbers either way, so neither form loses anything.

For this article, the important part is that a function and its list of wave amplitudes (its spectrum) hold the same information. The samples say what the value is at each place. The spectrum says how much of each wave the function contains. Either one determines the other, and the Fourier transform converts between them. It's worth converting because some problems are much easier on the spectrum side, and Poisson's equation is one of them.

Step 3 · Frequency spacePoisson's equation becomes division

Take a single wave of density, ρ = sin(kx). Differentiate it twice and you get −k² sin(kx), the same wave scaled by −k². A wave goes through the Laplacian with its shape unchanged, so for one wave the potential is that wave divided by −k². Any density is a sum of waves, and Poisson's equation is linear, so each wave can be solved on its own. The 262,144 coupled equations turn into 262,144 separate divisions.

Here it is one stage at a time: break the density into waves, divide, add them back up.

In code:

spectrum = fft2d(density);            // to frequency space
for (each wave k) {
  spectrum[k] /= -dot(k, k);          // this line solves Poisson's equation
}
potential = inverse_fft2d(spectrum);  // back to the mesh

In position space, Poisson's equation ties every cell to every other cell. In frequency space it's one division per cell. Dividing by k² also shrinks short waves much more than long ones, which is why the potential comes out so much smoother than the density. Periodic boundaries come for free, because the waves the FFT uses already wrap around the box. (The real shader also multiplies by a gentle Gaussian in k, which smooths out noise at the scale of single cells.)

Here it is running live. Drag to paint mass (shift-drag or right-drag erases). Every frame the demo solves Poisson's equation over its whole 256×256 mesh and draws the potential as a contour map, and 50,000 weightless tracer grains fall through it so you can see the slopes. The tracers lose a little speed every step so they settle into the wells instead of orbiting forever. Paint two blobs and their wells merge into one basin before the blobs touch. Paint near an edge and the well wraps around to the other side.

This is the third time in the series that a hard problem got easier by changing how it's represented. Part one grouped bodies into a tree so distant ones could be summarized. Part three sorted particles so neighbours sat next to each other in memory. Here the Fourier transform rewrites the density as waves, and in that form each wave can be solved by itself. Mathematicians call that diagonalizing the operator.

Step 4 · The FFTFrom n² to n log n

That leaves fft2d itself. Each output of a Fourier transform is a weighted sum over all n inputs, so computing all n outputs directly costs n², the same wall as in earlier parts. Straight from the definition, one 2D transform of a 512×512 mesh is 262,144² multiply-adds, about 69 billion. Each step needs two transforms (forward and back), and the cosmology demo takes two steps per frame.

The fast Fourier transform gets around this. A transform of size n splits into two transforms of size n/2, one over the even-numbered samples and one over the odd ones, recombined with one complex multiply per output. Do that recursively and the cost drops to n log n. The diagram below is the butterfly network for 16 points: log₂16 = 4 layers, where each layer pairs every value with a partner 8 positions away, then 4, then 2, then 1. It's the same layered shape as part three's prefix sum, with complex rotations in place of additions, and it maps onto a GPU just as well.

Below the diagram, two bars race a direct n² transform against an FFT on your CPU. Past 2,048 points the direct version would take too long, so it's timed at 2,048 and scaled up.

Step 5 · Particle-meshFrom particles to a field and back

The solver works on fields, but the matter is particles. The method that connects them is called particle-mesh, and each step goes like this:

deposit:   particles spread their mass onto the mesh   // particles → field
solve:     FFT → divide by −k² → inverse FFT            // potential
gradient:  difference neighbouring cells                // force at each cell
gather:    particles read the force where they are      // field → particles
kick and drift, wrap through the walls, repeat

Deposit and gather are the only parts with real design choices. The crude deposit drops each particle's whole mass into the cell it's in. That jitters: when a particle crosses a cell border its mass jumps from one cell to the next, the field jumps with it, and the particle feels the jump. The standard fix is cloud-in-cell (CIC). Each particle splits its mass over the four nearest cells, weighted by how much a cell-sized square centred on the particle overlaps each one, so the deposit moves as smoothly as the particle does.

The gather reads the force back with the same four weights. Depositing with one scheme and gathering with another would leave each particle feeling a small force from itself, and momentum wouldn't be conserved. Drag the particles below and compare CIC with the nearest-cell version:

That's the whole method: deposit, transform, divide, transform back, gather. No tree walk and no sort. Below is a cold blob like part four's ghost cloud, this time in a periodic box: 65,536 particles starting at rest, restarting every 14 seconds. The slider changes nothing about the particles, only the mesh the forces are computed on, from 64×64 to 512×512. On the coarse mesh the forces are blurry and the collapse stays smooth. On the fine one the same blob breaks up into much finer filaments and knots. The box is drawn tiled 3×3, so you can watch the neighbouring copies collapse in step, since they're the same particles.

That slider shows the method's main limit: below the size of one cell, there's no gravity to speak of. The mesh is the softening. For a galaxy that would be fatal. Part one needed forces down to the spacing between stars, which is why it used a tree. For cosmology it's a reasonable trade, because the structures of interest are large and fill the whole box. The tree spends bookkeeping to put resolution where matter clumps. The mesh gives every part of the box the same resolution and makes the whole box cheap.

Step 6 · Initial conditionsSeeds and an expanding box

The solver is done. The cosmology demo needs two more things.

The seeds. The young universe was almost featureless. The cosmic microwave background, mapped by COBE, WMAP and Planck, shows it was smooth to about one part in 100,000. The ripples are random, but they have a measured statistical fingerprint: how much ripple there is at each scale, called the power spectrum. So making initial conditions is short. Fill a grid with random noise, transform it, and scale each wave to match the spectrum you want. Then, instead of scattering particles randomly (which would add noise of its own), start them on a perfect lattice and displace each one along the resulting field, with a small velocity in the same direction. This runs on the CPU each time you re-seed, using three more FFTs.

The expansion. Space expands, and gravity has to work against it. Simulations fold the expansion into the coordinates. Comoving coordinates stretch along with space, so the box stays size 1 while real distances grow with the scale factor a(t). The expansion only shows up as powers of a in the kick and drift. Everything else stays the same: same kernels, same four FFT passes.

The clock runs in scale factor, from a = 0.02 (the universe at 2% of its present size) to a = 1, which is today. The ripples at the start are bigger than the microwave background's one part in 100,000, since they've had time to grow by then, but they're still small. The demo steps in equal increments of ln a, about 1,000 steps in all, so at 60 frames a second it reaches today in roughly eight seconds. Then the clock stops, but the particles keep moving and you can keep stirring. It uses a matter-only universe (called Einstein–de Sitter), which is simpler than the real one with its dark energy.

Run it below. The expansion on/off button shows why expansion matters. With it off, gravity runs unopposed and structure grows much faster, ending up as a few big clumps. With it on, gravity wins slowly and only where the seeds favoured it: knots first, then filaments between them, then empty voids as matter drains onto the web.

The ripple strength slider plays the part cosmologists give σ₈, the overall amplitude of the ripples: turn it down and structure shows up later and thinner. The tilt n slider sets how the ripples are split between large and small scales (the spectrum goes as kⁿ). Negative values give a few big filaments, and positive values give a foam of small clumps. Planck measured both numbers for the real universe, σ₈ at about 0.81 and the tilt of the primordial spectrum at about 0.965. They're two of the numbers that describe the standard cosmological model, ΛCDM, which has six free parameters in all. The sliders here are simplified stand-ins for them. The particle slider goes from 65,536 to just over a million.

One caveat: this box is 2D. Poisson's equation in 2D gives a force between two points that falls off as 1/r instead of 1/r². The kick and drift are set up so small ripples grow at the same rate as they would in a matter-only 3D universe, but once things get lumpy the details differ from the real thing. The method is the real one, though. In 3D, the mesh-based codes of the mid-1980s used the same ingredients: a mesh, FFTs and comoving kicks, plus a direct short-range correction (that version is called P³M). With 32,768 particles, Davis, Efstathiou, Frenk and White used them to make the case that cold dark matter, starting from ripples like these, could produce the clustering of galaxies we see.

EpilogueThree ways around n²

Across five parts the n² wall came down three times, with three different structures:

• The tree summarizes space. Far-away things blur together, so resolution can follow interest. Part one's pyramid answers "what's the total pull on me?" in log n steps, and pays in memory for the hierarchy.
• The sort organizes memory. Things that touch should sit next to each other in memory, not only in space. Part three's counting sort answers "who's next to me?" in constant time per particle, and pays with a scan.
• The transform changes the representation. Some equations are only hard in the form they were asked in. The FFT moves Poisson's equation to where it's a division, and pays n log n.

The big production codes use two at once. TreePM, the method in GADGET-2 that ran the Millennium Simulation's ten billion particles, splits gravity at a scale of a few mesh cells. The mesh handles the far field, periodic and cheap, and the tree handles the near field, sharp where matter clumps. Part four paired the tree with the sort. TreePM pairs the tree with the mesh.

Things I'd try next: a third dimension; finer meshes nested inside the deepest wells (adaptive mesh refinement); and initial conditions built from a real cosmological power spectrum instead of a tilt slider.

Next: Part six, The Fluid Instrument. The same Poisson solve, turned to pressure: a 3D fluid you can stir with your cursor or your hands.