Interactive 3D N-Body Gravity Simulation with Barnes-Hut Optimization

Why I Built This

I’ve always been fascinated by gravity. Not just the everyday kind that keeps us grounded, but the cosmic choreography that shapes entire galaxies. How do billions of stars dance together to form those beautiful spiral arms we see in Hubble images? What happens when two galaxies collide? These questions led me down a rabbit hole that eventually resulted in this interactive 3D simulation.

The N-body problem is basically this: if you have N objects all pulling on each other with gravity, how do they move? For two objects (like Earth and Moon), we can solve this exactly. But add just one more object and suddenly the math becomes… well, let’s just say there’s no clean solution. You need to simulate it step by step.

I wanted to show you two different approaches:

• Barnes-Hut Algorithm: A clever trick that makes the simulation blazing fast
• Direct Summation: The brute-force approach where every particle feels every other particle

But here’s what really excites me about this project…

The Real Story of Galaxy Formation

Before we dive into the simulation, let me clear up a common misconception. Galaxies didn’t form from explosions—that’s actually way cooler than that.

Picture this: About 13.8 billion years ago, the Big Bang created a universe that was almost perfectly smooth, but not quite. There were tiny density fluctuations—some regions had just a little bit more matter than others. These tiny differences were like seeds.

Over millions of years, gravity did its thing. The denser regions attracted more matter, growing slowly but steadily. Dark matter (which we can’t see but makes up most of the universe) clumped together first, forming invisible scaffolding called “dark matter halos.”

Then regular matter—hydrogen and helium gas—fell into these dark matter halos like water flowing into a bowl. As the gas collapsed, it started spinning (conservation of angular momentum), forming rotating discs. When the gas became dense and hot enough, stars began to form.

This is exactly what you’ll see in the simulation! Those rotating discs that emerge aren’t just pretty patterns—they’re recreating billions of years of cosmic evolution in fast-forward. The spiral arms that develop naturally from gravitational instabilities? That’s how we think our own Milky Way got its iconic shape.

The galaxy merger scenarios I included show another crucial part of the story. Galaxies aren’t isolated islands—they’re constantly interacting, merging, and growing. Most large galaxies, including ours, are thought to be the result of many smaller galaxies merging over cosmic time.

So no explosions—just gravity, time, and the universe’s incredible ability to build structure from simplicity.

This Is How We Study the Universe

When I started building this simulation, I didn’t realize I was essentially recreating the same methods that astronomers use to understand cosmic evolution. It turns out the algorithms running in your browser right now are the same ones powering some of the most ambitious scientific endeavors in human history.

Galaxy Formation in Your Browser

That rotating disc you’ll see in the simulation? That’s not just a pretty pattern—it’s how we think the Milky Way formed billions of years ago. Astrophysicists run massive versions of this simulation with millions of particles to understand how spiral arms emerge and evolve.

The “Galaxy Merger” preset I included recreates one of the most dramatic events in cosmic history. In about 4.5 billion years, our Milky Way will collide with Andromeda—and this is exactly how scientists study what will happen.

Cosmic Web Formation

The same Barnes-Hut algorithm I implemented here scales up to track dark matter particles forming the largest structures in the universe—the cosmic web of filaments that spans hundreds of millions of light-years. When you watch particles clumping together in the simulation, you’re seeing the same physics that created every galaxy cluster we observe.

A 1986 Revolution

Josh Barnes and Piet Hut published their breakthrough algorithm in 1986, and it completely transformed computational astrophysics. Before their work, simulations were limited to a few thousand particles—nowhere near enough to model realistic galaxies. Their clever insight about treating distant groups of particles as single masses suddenly made simulations with millions of particles possible.

Modern Supercomputing

Today’s astrophysics codes like GADGET, AREPO, and RAMSES use GPU acceleration to simulate systems with billions of particles. The Millennium Simulation—one of the most famous—followed 21.6 billion particles over the entire history of the universe, from shortly after the Big Bang to today.

What You’re Really Experiencing

When you play with this simulation, you’re not just watching pretty dots move around. You’re experiencing the same fundamental physics that governs:

• How the cosmic web formed after the Big Bang
• Why our galaxy has spiral arms
• What will happen when Milky Way and Andromeda collide
• How globular clusters evolve over billions of years

Every time you adjust the parameters and watch the system evolve, you’re doing the same kind of experiments that help us understand our place in the cosmos.

The Physics Behind the Simulation

Gravitational N-Body Dynamics

Each particle in our simulation follows Newton’s laws under gravitational forces from all other particles. The gravitational force between two particles with masses $m_i$ and $m_j$ separated by distance $r_{ij}$ is:

Where:

• $G$ is the gravitational constant
• $\vec{r}_{ij} = \vec{r}_j - \vec{r}_i$ is the separation vector
• The negative sign indicates attraction

The total force on particle $i$ is the sum over all other particles:

This leads to the equation of motion:

The Computational Challenge

For N particles, computing all pairwise forces requires $\frac{N(N-1)}{2} \approx \frac{N^2}{2}$ calculations per time step. This O(N²) scaling becomes prohibitively expensive for large systems:

• 1,000 particles: ~500,000 force calculations
• 10,000 particles: ~50,000,000 force calculations
• 100,000 particles: ~5,000,000,000 force calculations

Real astrophysical simulations often involve millions or billions of particles, making direct methods impractical.

The Barnes-Hut Algorithm

Core Principle

The Barnes-Hut algorithm, developed by Josh Barnes and Piet Hut in 1986, reduces computational complexity from O(N²) to O(N log N) using a clever approximation: distant groups of particles can be treated as single massive particles located at their center of mass.

Spatial Decomposition with Octrees

The algorithm works by:

1. Spatial Subdivision: Recursively divide 3D space into octants (8 cubic regions)
2. Hierarchical Tree: Build an octree where each node represents a spatial region
3. Mass Distribution: Store total mass and center of mass for each node
4. Multipole Approximation: Use the multipole acceptance criterion to decide when to approximate

The Multipole Acceptance Criterion

For each particle, we traverse the octree and at each node ask: “Can I treat this entire subtree as a single particle?”

The criterion is: $\frac{s}{d} < \theta$

Where:

• $s$ = size of the tree node (spatial extent)
• $d$ = distance from particle to node’s center of mass
• $\theta$ = theta parameter (accuracy threshold)

If the criterion is satisfied: Treat the entire node as a single particle If not: Descend into child nodes and repeat

Theta Parameter Trade-off

• Small θ (0.3-0.5): High accuracy, more computations, approaches direct method
• Large θ (1.0-1.5): Fast computation, lower accuracy, more approximation error

Typically θ ≈ 0.5-0.8 provides good balance between speed and accuracy.

Numerical Integration: Leapfrog Method

We use the leapfrog integration scheme, which is symplectic (energy-conserving) and second-order accurate:

1. Kick (half step):    v(t+dt/2) = v(t) + a(t) × dt/2
2. Drift (full step):   x(t+dt) = x(t) + v(t+dt/2) × dt
3. Compute forces:      a(t+dt) = F(x(t+dt))/m
4. Kick (half step):    v(t+dt) = v(t+dt/2) + a(t+dt) × dt/2

This method preserves energy much better than simple Euler integration, crucial for long-term stability.

Simulation Setup: Rotating Galactic Disc

Our initial conditions create a rotating disc that mimics galaxy formation:

Position Distribution

Particles are distributed with radius $r \propto \sqrt{\text{random}}$ to create realistic density profiles, with small vertical dispersion to form a thin disc.

Velocity Profile

Each particle receives approximately circular velocity plus perturbations:

• Circular velocity: $v_{\text{circ}} = \sqrt{\frac{GM_{\text{enc}}(r)}{r}}$
• Bar-like shear: Creates spiral arm instabilities
• Random velocities: Adds realistic velocity dispersion

This setup naturally evolves into spiral galaxy structures through gravitational interactions.

Parameter Guide

Theta (θ): Barnes-Hut Accuracy

• Range: 0.3 - 1.5
• Effect: Controls speed vs. accuracy trade-off
• Low values (0.3-0.5): More accurate, slower computation
• High values (1.0-1.5): Faster computation, more approximation
• Recommended: 0.5-0.8 for good balance

Time Step (dt): Integration Precision

• Range: 0.001 - 0.03
• Effect: Smaller steps = better accuracy but slower simulation
• Too large: Energy drift, numerical instabilities
• Too small: Unnecessarily slow computation
• Recommended: 0.01 for stable evolution

Softening (ε): Gravitational Smoothing

• Range: 0.002 - 0.05
• Purpose: Prevents singular forces when particles get very close
• Physics: Represents finite particle size or quantum effects
• Effect: Larger values make forces smoother but less realistic
• Recommended: 0.01 for stable dynamics

Particles (N): System Size

• Range: 1,000 - 8,000
• Performance:
  • Barnes-Hut scales as O(N log N)
  • Direct method scales as O(N²)
• Memory: Higher N requires more RAM
• Visual: More particles show finer structure

What to Observe

Spiral Arm Formation

Watch how the initially smooth disc develops spiral density waves through gravitational instabilities. This mirrors real galaxy evolution.

Performance Comparison

Switch between Barnes-Hut and Direct methods to see the dramatic performance difference. With 3000+ particles, Barnes-Hut runs ~10-100× faster while maintaining good accuracy.

Energy Conservation

Monitor the energy drift graph. Good simulations maintain energy to within ~1-5%. Larger drift indicates numerical problems or inadequate time resolution.

Total energy: $E = K + U$

Kinetic energy: \(K = \frac{1}{2}\sum_{i=1}^{N} m_i |\vec{v}_i|^2\)

Potential energy: \(U = -\frac{1}{2}\sum_{i=1}^{N}\sum_{j \neq i} \frac{G m_i m_j}{|\vec{r}_{ij}|}\)

Energy drift: \(\Delta E = \frac{E(t) - E(0)}{|E(0)|} \times 100\%\)

Good simulations maintain \(\|\Delta E\| < 1\) over hundreds of dynamical times.

Dynamic Range

Observe both:

• Large-scale structure: Overall spiral pattern evolution
• Small-scale dynamics: Individual particle orbits and interactions

Ready to Play God with Gravity?

I’ve spent way too many hours tweaking this simulation (my browser history can confirm), and now it’s your turn! The controls below let you experiment with the same physics that shapes galaxies.

Here’s what I love doing:

• Start with the “Single Disc” and watch spiral arms emerge naturally—it’s mesmerizing
• Switch to “Galaxy Merger” and witness cosmic destruction in real-time
• Crank up the particle count and watch Barnes-Hut work its magic (try comparing it with Direct mode… if you dare)
• Play with the theta parameter—you’ll see the accuracy vs. speed trade-off in action

The energy drift graph at the bottom is my favorite debugging tool. If it stays flat, you’re doing physics right. If it goes crazy, well… that’s when things get interesting.

Go ahead, break the universe. I’ve already done it countless times while building this!