kelvin@cluxnei:~/blog$ cat simulating-the-universe.md

2026-09-12 EN ~25 min read ler em português

Simulating the universe: the physics behind 700 bodies in a browser tab

This is an exploration, not a story. We are going to open up a gravitational simulation and look at every part: the equations and what each one is for, the real names of the techniques, the measured performance numbers, and the part almost nobody writes down, which is where the simulation lies on purpose, with the cost of each lie calculated.

The physics, from the start

Before looking at any code, it is worth assembling the whole vocabulary. There are few formulas and each one answers a specific question.

The only force here

Everything starts with Newton's law of universal gravitation. Two bodies of mass m₁ and m₂, separated by a distance r, attract each other with force:

Two facts turn that simple formula into a hard problem. The first is that the force falls with the square of the distance but never reaches zero: every body pulls every body, forever. The second is that acceleration is force divided by mass, so the acceleration a body i feels is the sum of the contributions from all the others:

With two bodies this has an exact solution, which is Kepler's ellipses. With three or more it does not. Poincaré showed in the late nineteenth century that no general closed form solution exists for the N body problem. That is why you simulate: not because the mathematics is unknown, but because the only way to know where the bodies will be a year from now is to advance time in small steps and add.

Density, radius and the derived quantities

A body does not get an arbitrary radius. Given mass and density, the radius follows from assuming a sphere:

Tying radius to density is the difference between a model and an animation: from here on, collisions, tides and densities start meaning something.

The temperature at the centre of a star comes out of an equilibrium argument: the gas pressure has to hold up the weight of the layers above. From that comes an estimate known as the virial temperature:

Applied to the Sun it gives 1.39 × 10⁷ K against the real 1.57 × 10⁷ K, so 12% error out of a one line calculation. That is good enough to decide what burns and what does not.

The distances that decide the behaviour

Three radii show up constantly, and each one marks the boundary of a different regime.

The Hill sphere is the region where a body's own gravity wins against the star it orbits. It is its territory:

The tidal radius is where the difference in gravity between the near and far sides of a body beats the gravity holding it together. Past that, it comes apart:

And the Schwarzschild radius is a black hole's event horizon, the distance at which escape velocity equals the speed of light:

Hold on to that last one. Comparing r_s against the simulation's other scales is what will decide, later, what is worth simulating and what is decoration.

Units: why G = 4π²

In SI, the gravitational constant is 6.674 × 10⁻¹¹ m³·kg⁻¹·s⁻². It is an awful number to work with: masses land around 10³⁰, distances around 10¹¹, and a float64 spends a good part of its precision just carrying exponent.

The classic way out in astronomy is to change units. Measure mass in solar masses, distance in astronomical units and time in years, and something lovely happens.

Consider a circular orbit. The centripetal acceleration has to be exactly gravity:

Now impose the definition of the astronomical unit itself: a body at 1 AU from a 1 M☉ star completes its orbit in 1 year. That orbit's speed is the circumference over the period, which is 2π. Substituting r = 1 and M = 1:

And that hands you, for free, the best test this simulation has. It is written as doctrine at the top of units.js:

Put one body at 1 AU with speed 2π AU/yr and it must close its orbit in exactly one year.

It is a falsifiable test that does not depend on any other part of the code being correct. If it fails, the foundation is wrong and nothing else is worth looking at.

The singularity and Plummer softening

Newton's law has a fatal defect for anyone simulating it: as r → 0, the force goes to infinity. Two bodies passing very close get an absurd acceleration in a single timestep and shoot off, and the simulation explodes. That is not an implementation bug. It is the formula being honest about a situation the timestep cannot resolve.

The naive fix is to clamp acceleration at an arbitrary ceiling. It works in the sense that nothing flies off to infinity, but it is physically indefensible: every time the clamp fires you are injecting or removing energy from the system without anything physical having happened.

The fix used here is Plummer softening, which adds a length ε to the denominator:

There is no discontinuity, no energy created from nothing, and it even saves a vector normalisation that could produce NaN when the distance is exactly zero.

The criterion for choosing ε is easy to state: it has to be far smaller than the orbital separations you care about, or it would distort the very orbits you want to measure. A fifth of a solar radius against separations of order 1 AU clears that comfortably. The measured effect on a test orbit closing is negligible: phase error 1.030 × 10⁻³ rad with softening on against 1.011 × 10⁻³ rad with it off.

The integrator, and how to know it is right

Integrating means: given where the bodies are and the acceleration they feel, work out where they will be an instant later. The simplest method is Euler, which advances position using current velocity and velocity using current acceleration. It is first order, and in practice orbits visibly precess and energy drifts.

The method used here is Velocity-Verlet, which costs almost the same, since you only have to keep the previous acceleration, and is second order:

Beyond the order, it has a property that matters a lot in gravitation: it is symplectic. Energy does not drift systematically up or down, it oscillates around the right value. A closed orbit stays closed after thousands of laps.

The part that actually proves something

"It looks like it works" is not verification. The right way to test an integrator is to measure its order of convergence: if it is second order, dividing the step by k has to divide the error by k².

step reduced by:   35×
error reduced by:  1222×
expected (35²):    1223

The integrator has exactly the order of convergence theory predicts. That is a claim which could have come out false, and that is what makes it useful.

The timestep belongs to the scenario, not to the program

An engineering detail with a direct physical consequence: there is no global timestep. Every scenario publishes its own, because a tight binary needs a far smaller one than a 20 AU disc, and reusing the disc's would visibly break the binary.

The criterion comes from the innermost orbit, which is always the fastest. In the default disc it takes about 0.35 years, and Velocity-Verlet wants well over 20 steps per orbit:

FIXED_DT = 0.0035 yr   about 100 steps on the inner orbit

It is also why the black hole scenarios run so much slower in your browser: near a 10⁶ M☉ body the required step is tiny, and each second of wall clock covers far less simulated time.

Barnes-Hut: trading accuracy for time

Computing the force between every pair of n bodies costs O(n²). With 800 bodies that is 319,600 pairs per step, and at 100 steps per orbit it adds up.

The Barnes-Hut algorithm fixes this with a simple observation: a distant cluster of bodies pulls about the same as a single point mass at its centre of mass. Space is divided into an octree, where each cube subdivides into eight until one body is left per cell, and when computing the force on a body you walk the tree deciding cell by cell with the opening criterion:

That drops the cost to O(n log n). The price is error, and the error is measured, not estimated:

That last row is another falsifiable test: with θ going to zero, the tree stops approximating anything and has to reproduce the direct sum. It does, to fifteen digits.

The uncomfortable honesty

Barnes-Hut is the most cited optimisation in N body simulation. At this scale, though, the measured gain was:

800 bodies:    ~1.1×   the crossover sits right about here
~2000 bodies:  ~2.2×   from here it genuinely pays

At the default body count, the tree barely pays for its own construction. It is there because the body count is adjustable and the user can ask for 4000, but it would be dishonest to sell it as the reason the simulation is fast.

Why every star is summed exactly

Here is a decision worth more than the whole optimisation. The Barnes-Hut tree is applied only to the light bodies. Every star is summed exactly, body by body, outside the tree.

The reasoning: the stars carry almost all the mass of the system. Each light body's orbit is essentially its orbit around those few terms. Approximating the stars with a multipole would put the error exactly where it dominates the result, which is saving in the wrong place. And since there are only a handful to a few hundred stars, the cost is O(S·n) and exact, against O(n log n) and approximate for everything else.

The two passes touch disjoint pairs: the tree never sees a star, and the direct sum covers every pair involving one. Nothing is counted twice, nothing is missed.

A black hole is on that list too, and the comment in the code makes a point of saying this is not a claim that it shines. It is simply the direct summation set, and a hole is the heaviest object in any system containing it.

The momentum the tree breaks, and the fix

There is a subtle side effect. The forces computed by the tree are not pairwise symmetric: body A feels the cell containing B, but B feels a different cell containing A. Newton's third law stops holding exactly, and a spurious net force is left on the whole system.

The fix is elegant. Sum that net force and subtract it as a uniform acceleration from every body that went through the tree. This restores the condition below without altering any relative acceleration, and therefore without touching any orbit:

raw net force from the tree:   6 × 10⁻³
after the correction:          3.9 × 10⁻¹⁶

measured with 800 bodies, after 531 mergers

Composition: alpha chain fusion

Every body carries a chemical composition. Not a label, but mass fractions of eight species of the alpha chain, the fusion sequence that governs a star's life:

H → He → C, O → Ne, Mg → Si → Fe

each step needs a higher central temperature than the one
before. Iron is the end of the line, because fusing iron
consumes energy instead of releasing it.

Each body's central temperature is recomputed every step with the virial formula from the introduction, and it decides which step is lit. This is an equilibrium, not an integration: a star's thermal time and its dynamical time differ by orders of magnitude. The Sun's Kelvin-Helmholtz time is tens of millions of years, its dynamical time is half an hour. At the scale of the step, the thermal structure is quasi-static, and computing the equilibrium is cheaper and more stable than evolving it.

The bug only a test would catch

The chain as first specified had absorbing states: nitrogen and magnesium had no way out. A core at 5 × 10⁹ K, hot enough to reach iron, finished with 69% magnesium and stopped there.

The fix was adding the two missing transitions: N → O, since nitrogen is a catalyst of the CNO cycle but is destroyed by alpha capture once helium ignites, and Mg → Si. The comment left in the code sums up the invariant:

Iron is the only absorbing state in this chain, which is the point.

Another case of the same kind: the ignition threshold, if it were temperature alone, would light brown dwarfs. The virial estimate runs about 4× high near the limit, because in very small objects what holds the star up is electron degeneracy pressure, which costs no heat. Ignition is now conditioned on mass, and the result is correct:

0.079 M☉  →  brown dwarf, T_c = 1.22 × 10⁷ K, burns nothing
0.080 M☉  →  star, burns

The poetic licences

This is the section I most wanted to write, because it is the one usually left out. Every simulation lies somewhere. The difference between a model and an animation is whether the lies are marked and whether their cost was calculated.

1. The capture radius, and the project's largest lie

At true scale, a planet has a radius of order 10⁻⁵ AU against orbital separations of order 1 AU. The geometric cross section is so absurdly small that a real disc takes 10⁷ to 10⁸ years to build planets. In a browser tab, nothing would ever happen.

So the capture radius is multiplied by 800. That compresses planet formation into a few thousand simulated years. What saves this from being a trick is what comes with it:

• The dynamics is not touched. Orbits, resonances and scattering all stay exact. Only the target size is exaggerated.
• The cost was measured: at factor 800, an embryo's capture radius at 1 AU is 3.3 times its own Hill sphere. That is, inner disc bodies capture each other well outside the territory where they should rule.
• The honest alternative is written right beside it: at factor 250 the capture radius stays inside the Hill sphere throughout the inner disc, at the cost of 11 Earth masses in 150 years instead of a super-Earth in 27.

And the boundary between licence and physics is marked in one line of the code. Gravitational focusing, which enhances the cross section by a factor (1 + v_esc²/v_rel²) and matters when the relative velocity is low, is real physics, not licence.

There is also a lie that was consciously not committed: the central star absorbs at 10 solar radii, not 800. The comment explains why, since 800× would be 3.7 AU and it would eat the whole disc.

2. The drawn size

One solar radius is 0.00465 AU. In a 20 AU disc, a true-scale planet takes up a few millionths of the scene: the screen would be black. Bodies are drawn inflated, on a compressive curve with exponent 0.55, so that a Jupiter does not dwarf an Earth by the factor it really would.

That scale was tuned by measurement, not by taste. With the value that seemed reasonable at first, the entire dynamic range collapsed:

Tidal disruption: the test against nature

When a star passes too close to a supermassive black hole, the difference in gravity between its near and far sides beats its own self gravity and it comes apart. Half the debris stays bound and returns, the other half escapes. It is a real event, observed dozens of times, and it has a characteristic signature.

The design question was whether this could be simulated without hydrodynamics. The answer came from comparing the tidal radius with the softening length, that ruler of honesty from the earlier section:

At a supermassive hole, the region where the star comes apart is hundreds of times larger than the scale at which the simulation already stopped resolving. Which means there is no need to fake it. On crossing the tidal radius, the star is replaced by 600 bodies carrying its mass and composition, and the physics that already existed does the rest. The leading particles orbit faster than the trailing ones, and the stretching appears on its own.

Spaghettification stops being drawn and becomes a consequence.

It did not work at first, and the reason matters

On the first attempt the fragments were laid out along a solar radius. That put them 1.5 × 10⁻⁵ AU apart, 65 times closer than the softening length. And below the softening, the Plummer kernel is not gravity: it is a linear spring.

The result was immediate and unambiguous. The energy distribution collapsed into four spikes in under 100 steps, and the fitted exponent came out −0.41, with r² = 0.04. That was not an inaccuracy. It was the model describing something else.

The fix was widening the stream until the mean spacing was a few softening lengths. And it was classified, in the code itself, as poetic licence of the same species as the 800 factor, because that is exactly what it is.

The result

The observational signature of a tidal disruption event is the debris return rate following a power law:

The exponent comes out about 2.5% steeper than theory, consistently, and the cause was diagnosed rather than swept away: late returns are missing at the end of the window, which tilts the fit. A purely Keplerian fit over the initial energies gives −1.663.

Where it does not work

With a 10 M☉ stellar mass black hole, conservation stays exact and the split is still half and half, but the physical result is no good. At that mass, the star's own absorption radius is five times larger than the tidal radius, so the "disruption" fires at nine tidal radii, nowhere near where it should. The useful crossover sits around 10⁴ M☉.

What was deliberately refused

The most mature decision in the whole project is something that was not implemented.

Near a black hole, Newtonian gravity is wrong. There is a well known shortcut for capturing the relativistic effects without solving Einstein, the Paczyński-Wiita pseudo-Newtonian potential:

It would be easy to add and would make a good line in a README. It was refused, and the reason is recorded in the code with an explicit request that nobody add it later by accident. The arithmetic:

10 M☉ hole:
  r_s          = 1.97 × 10⁻⁸ AU
  ISCO (3·r_s) = 6 × 10⁻⁷ AU
  softening    = 1 × 10⁻³ AU

the ISCO sits 4 to 5 orders of magnitude BELOW the scale
at which the force law is already deliberately wrong

It would be decoration painted over a place the simulation does not resolve, at the cost of a special case in the hottest loop in the program and of exact momentum conservation. For comparison: at a 10⁶ M☉ hole the r_s / softening ratio is 19.7, and there the horizon is resolvable. That is why supermassive black holes got in and relativity did not.

By the same criterion, real hydrodynamics (SPH) stayed out. It would need pressure, shocks and an equation of state, and would change the character of the whole project. The disruption fragments are gravitational masses, and the scenario's own text says so to the user.

The summary

If there is one thing to take from this text, let it be the criterion that shows up in every section: choose the scale at which you are willing to be wrong, and then refuse everything that happens below it.

The softening length defines that scale. It justifies the softening, it justifies refusing relativity, it condemns the first version of the debris stream and it decides which black holes are worth having. A single number, applied consistently, does more for a model's honesty than any amount of physics piled on top.

The simulation lives at cluxnei.dev/universe-simulation/ and runs entirely in your browser, no server and no install. The source is on GitHub, including every comment quoted here, which are the best part.

kelvin@cluxnei:~/blog$