Galactic N-Body Simulation

Computational Astrophysics

Newtonian Gravity

The core of the simulation integrates Newton's law of universal gravitation across an N-body system. To prevent infinite forces and numerical instability during close encounters (singularities), a softening parameter $\epsilon$ is introduced.

The acceleration $\mathbf{a}_i$ of particle $i$ due to all other particles $j$ is calculated as:

$$ \mathbf{a}_i = G \sum_{j \neq i} \frac{m_j}{(\|\mathbf{r}_{ij}\|^2 + \epsilon^2)^{3/2}} \mathbf{r}_{ij} $$

Where $\mathbf{r}_{ij} = \mathbf{r}_j - \mathbf{r}_i$ is the distance vector, $m_j$ is the mass of the interacting particle, and $G$ is the gravitational constant.

Isothermal Dark Matter Halo

Visible matter alone cannot account for the observed flat rotation curves of galaxies. To model this without the computational overhead of millions of unseen particles, the simulation implements an implicit Isothermal Dark Matter Halo force.

This adds an inward acceleration component directed towards the galactic center, parameterized by:

$$ a_{dm} = -\frac{V_{flat}^2 \cdot r}{r^2 + \epsilon_{core}^2} $$

Where $V_{flat}$ represents the asymptotic circular velocity governed by the dark matter strength, and $\epsilon_{core}$ is the core radius of the halo.

Salpeter Initial Mass Function (IMF)

The mass distribution of stars in the simulation follows the realistic Salpeter Initial Mass Function, rather than a uniform or arbitrary assignment. The IMF describes the empirical distribution of initial masses for a population of stars.

$$ \xi(M) = \xi_0 M^{-2.35} $$

This power-law distribution means low-mass stars are exponentially more common than massive stars. In the simulation, this continuous mass spectrum is dynamically sorted to differentiate between heavy "active" stars (which exert gravity) and light "passive" stars (which only receive gravity).

Stellar Color & The H-R Diagram

Star colors in the simulation are not randomly assigned; they are correlated with a star's mass, mapped roughly to the main sequence of the Hertzsprung-Russell (H-R) diagram.

More massive stars burn hotter and brighter, shifting their emission spectra towards the blue end, while less massive stars are cooler and emit primarily in the red spectrum. The simulation visually distinguishes stars by rendering the most massive active bodies with intense blue-white hues, grading down to yellow, orange, and deep red for the numerous lightweight passive stars.

Numerical Analysis

Symplectic Leapfrog Integration

Standard explicit Euler integration diverges rapidly for orbiting bodies, violating the conservation of energy. To solve this, the simulation utilizes a 2nd-order Symplectic Leapfrog Integrator. By staggering the velocity and position updates by a half-step $\Delta t / 2$, this method maintains long-term phase-space area and ensures orbital stability over thousands of revolutions.

$$ \mathbf{v}_{i}(t + \frac{\Delta t}{2}) = \mathbf{v}_{i}(t - \frac{\Delta t}{2}) + \mathbf{a}_{i}(t) \Delta t $$ $$ \mathbf{r}_{i}(t + \Delta t) = \mathbf{r}_{i}(t) + \mathbf{v}_{i}(t + \frac{\Delta t}{2}) \Delta t $$

Leapfrog integration guarantees time-reversibility and reliable energy conservation over long durations, making it a critical algorithm for high-fidelity galactic dynamics where performance and accuracy must be strictly balanced.

Algorithmic Architecture

Brute Force Gravity (O(N²))

The most accurate, yet computationally expensive, method for calculating gravitational forces is the Direct N-Body (Brute Force) approach. In this model, every particle calculates its gravitational interaction with every other particle in the system.

While this ensures perfect accuracy, the time complexity scales quadratically at $O(N^2)$. Simulating 10,000 particles requires 100,000,000 force calculations per frame, making it intractable for macroscopic galactic simulations on typical hardware, but useful as a baseline truth for smaller clusters.

Barnes-Hut & QuadTree (O(N log N))

To overcome quadratic scaling, the simulation implements the Barnes-Hut algorithm. This spatial partitioning method divides the 2D simulation space into a hierarchical QuadTree. Particles are recursively grouped into hierarchical nodes (cells).

If a group of distant particles (a node) is sufficiently far away from the particle being evaluated, the algorithm approximates the entire group as a single center of mass rather than calculating individual forces. This drastically reduces the time complexity to $O(N \log N)$, enabling the real-time simulation of significantly larger particle counts without a noticeable loss in macroscopic accuracy.

Active/Passive Computational Subsetting

Even with Barnes-Hut, rendering a visually dense galaxy of 100,000+ stars on CPU threads is challenging. To solve this, we introduce an Active/Passive heuristic, splitting the galaxy into two mass-distinct sets:

• Active Set ($M_{heavy}$): The top percentage of most massive stars. These stars are inserted into the QuadTree and exert gravity on all other stars.
• Passive Set ($M_{light}$): The vast majority of low-mass stars. These stars possess negligible mass compared to the galactic core. They do not exert gravity and are not inserted into the QuadTree, but they do receive gravitational pull from the active set.

This hybrid approach effectively reduces the computational complexity from $O(N \log N)$ to $O(N \cdot k)$ where $k \ll N$, allowing for massive visual density while preserving accurate macroscopic galactic dynamics.

High Performance Computing

CPU Physics Engines & Threading

The simulation provides two distinct CPU-bound physics engines, each optimized for different classical N-body algorithms.

• Brute Force Engine: Implements an $O(N^2)$ direct-sum calculation using Velocity-Verlet integration. It maximizes efficiency on a single thread by exploiting Newton's Third Law (calculating the force between $i$ and $j$ once and applying it to both), effectively halving the required calculations.
• Barnes-Hut Engine: Utilizes a QuadTree data structure to reduce the computational complexity to $O(N \log N)$. This engine dynamically reconstructs the spatial tree every frame and evaluates forces based on a distance threshold $\theta$.

To prevent the physics calculations from blocking the main UI thread, the engines are executed asynchronously. A WorkerBridge spawns a dedicated Web Worker (physics.worker.ts). Memory is not copied between threads; instead, a SharedArrayBuffer is utilized to provide zero-copy concurrent memory access. The main thread and the worker synchronize state atomically using Atomics.wait and Atomics.notify, ensuring strict deterministic execution without memory locking overhead.

GPU Physics & Compute Architecture

For maximum performance and particle density, the simulation delegates both physics integration and rendering to the Graphics Processing Unit via the WebGPU API.

• Compute Shaders: Gravitational forces and velocity updates are calculated in parallel across thousands of GPU cores using WGSL compute shaders (shaders.wgsl). The compute pipeline employs a "ping-pong" buffer architecture: in frame $A$, it reads state from Buffer 1 and writes the updated state to Buffer 2. In frame $B$, it reverses the process. This guarantees that all instances read from a consistent state without complex thread synchronization.
• Hardware Optimizations: The gravity kernel leverages specialized hardware instructions, such as inverseSqrt, to rapidly compute the inverse distance required for the gravitational force vector, significantly outperforming standard CPU math libraries.
• Unified Render Pipeline: Unlike typical implementations, the calculated particle positions and velocities are never read back to the CPU memory (preventing costly PCI-e bus transfers). Instead, the GPU's Render Pipeline directly reads the compute shader's output buffer (Storage Buffer) as Instance Data. Particles are drawn as heavily optimized instanced Quads, with perfect circular appearances achieved via fragment shader pixel discarding and mathematical gradients.