3D N-Body Systems:
Octree Dynamics & Collision Detection

In astrodynamics and scientific computation, the N-body problem requires solving the equations of motion for $N$ particles interacting under physical forces (like gravity or electrostatics). In a full three-dimensional space, the gravitational force vector $\mathbf{F}_i$ acting on a particle $i$ due to all other particles is formulated as:

Where $\mathbf{r}_{ij} = \mathbf{x}_j - \mathbf{x}_i$ is the displacement vector, $m_i$ and $m_j$ are the masses, $G$ is the gravitational constant, and $\epsilon$ is a softening factor inserted to prevent numerical singularities as the inter-particle distance approaches zero.

As the particle count $N$ scales, computing all pairwise forces directly has a computational complexity of $\mathcal{O}(N^2)$. Moving from 2D coordinates to 3D trajectories substantially increases the computational burden: evaluating square roots for distance calculations, handling vector projections on three coordinate axes, and computing complex spatial collisions require a structured approach to spatial partitioning.

To bypass the $\mathcal{O}(N^2)$ pairwise calculation bottleneck, we partition space hierarchically. In a three-dimensional domain, the direct equivalent of a 2D Quadtree is an Octree.

An Octree begins with a single cubical bounding box enclosing the entire simulation space. As particles are inserted, any node containing more than a threshold capacity (typically one particle) splits into eight child octants along its spatial midplanes. The eight children represent the combinations of:

• Left / Right (negative/positive $X$ axis)
• Bottom / Top (negative/positive $Y$ axis)
• Back / Front (negative/positive $Z$ axis)

This subdivision creates an 8-way branching search tree. When traversing the tree to calculate gravitational forces or detect collisions, we can ignore entire regions of space that are far away or do not overlap with our target coordinates.

While both trees are hierarchical spatial data structures, the transition from two to three dimensions shifts the branching factor and memory layout.

Although the asymptotic time complexities look similar, the constant factors in 3D are significantly higher. Constructing an Octree requires partitioning space across three planes ($xy$, $yz$, and $xz$), resulting in more expensive coordinate checks per node insertion. However, because the branching factor is 8, the depth of an Octree for a uniform distribution of particles is shallower than that of a Quadtree:

In a particle simulation, check-looping every particle against every other to see if they overlap results in a costly $\mathcal{O}(N^2)$ double loop. By utilizing the Octree, we can prune checks by querying only the nodes that overlap a target particle's radius.

We approximate this via Axis-Aligned Bounding Box (AABB) queries. For a particle located at $(x, y, z)$ with collision radius $R$, we construct a query box:

We traverse the Octree starting at the root. If the query box does not overlap a node's bounding box, we discard that node and all of its descendants immediately. If it does overlap, we recurse into its children. Once we reach a leaf node containing particles, we perform exact Euclidean distance checks:

Below is a clean python reference demonstrating how to split a 3D bounding box into eight child octant coordinates during Octree insertion:

class BoundingBox3D:
    def __init__(self, min_x, min_y, min_z, max_x, max_y, max_z):
        self.min_x, self.min_y, self.min_z = min_x, min_y, min_z
        self.max_x, self.max_y, self.max_z = max_x, max_y, max_z

    def split(self):
        """Splits current 3D box into 8 child bounding boxes (octants)"""
        mid_x = (self.min_x + self.max_x) / 2
        mid_y = (self.min_y + self.max_y) / 2
        mid_z = (self.min_z + self.max_z) / 2

        return [
            # Bottom-Left-Back to Top-Right-Front
            BoundingBox3D(self.min_x, self.min_y, self.min_z, mid_x, mid_y, mid_z), # 0
            BoundingBox3D(mid_x, self.min_y, self.min_z, self.max_x, mid_y, mid_z), # 1
            BoundingBox3D(self.min_x, mid_y, self.min_z, mid_x, self.max_y, mid_z), # 2
            BoundingBox3D(mid_x, mid_y, self.min_z, self.max_x, self.max_y, mid_z), # 3
            BoundingBox3D(self.min_x, self.min_y, mid_z, mid_x, mid_y, self.max_z), # 4
            BoundingBox3D(mid_x, self.min_y, mid_z, self.max_x, mid_y, self.max_z), # 5
            BoundingBox3D(self.min_x, mid_y, mid_z, mid_x, self.max_y, self.max_z), # 6
            BoundingBox3D(mid_x, mid_y, mid_z, self.max_x, self.max_y, self.max_z)  # 7
        ]