ArborX: Portable Geometric Search Library

• ArborX is a high-performance geometric search library built on Kokkos, designed for efficient range, k-nearest-neighbor, collision, and ray queries.
• It evolved from a narrow BVH wrapper to a flexible framework with a unified C++20 API, supporting callbacks, distributed search, and advanced algorithms like DBSCAN and EMST.
• The library achieves significant speedups through data-parallel BVH construction, stackless traversal, and hardware-specific optimizations across CPUs and GPUs.

ArborX is a performance portable geometric search library built on top of Kokkos, originally introduced as an open-source C++ library for modern supercomputing architectures and later extended to support exascale applications and a broader algorithmic surface in version 2.0. Its core abstraction is a bounding-volume-hierarchy (BVH) index over geometric primitives, used to accelerate range, $k$-nearest-neighbor, collision, and ray queries, while preserving portability across host and device back ends through Kokkos execution and memory spaces. In the 2.0 release, ArborX exposes a unified C++20 API that runs efficiently on NVIDIA CUDA, AMD HIP, Intel SYCL, and multicore CPUs with OpenMP or Pthreads, and it adds a generalized interface, callbacks, distributed search, brute-force fallback, ray tracing, clustering, Euclidean minimum spanning tree (EMST), and moving-least-squares (MLS) interpolation (Lebrun-Grandié et al., 2019, Prokopenko et al., 2024, Prokopenko et al., 31 Jul 2025).

1. Origins, objectives, and application domain

ArborX was developed to address a specific combination of requirements that earlier C++ spatial-index libraries did not target simultaneously: low-overhead parallel tree build, batched query execution on millions of points with predictable memory usage, and performance portability across multi-core CPUs, many-core GPUs, and future architectures (Lebrun-Grandié et al., 2019). The 2019 description places these requirements in scientific applications such as computational mechanics, multiphysics coupling, computer vision, cosmology, mesh searches, particle-in-cell, SPH, and Lagrangian-Eulerian couplers, where proximity search is a recurring kernel and brute-force all-pairs comparisons cost $O(n^2)$ (Lebrun-Grandié et al., 2019).

The library’s later evolution was strongly shaped by exascale cosmology. The 2024 exascale paper describes ArborX as developed as part of the Exascale Computing Project (ECP), and it ties major enhancements to a collaboration with the HACC cosmology code, where in-situ halo finding is implemented with DBSCAN and each MPI rank handles hundreds of millions of particles (Prokopenko et al., 2024). That paper emphasizes that the pressure to support exascale supercomputers from different vendors required performance portability not only for the search index itself but also for the clustering workflow layered on top of it (Prokopenko et al., 2024).

Version 2.0 broadens the scope beyond the original BVH-centric interface. The technical report characterizes the release by five design goals: generalized interface, callbacks, distributed search, expanded algorithms, and performance improvements (Prokopenko et al., 31 Jul 2025). A common misconception is that ArborX is only a BVH wrapper for 3D axis-aligned boxes. That description matches the 1.x API, but the 2.0 release explicitly expands the admissible value types, bounding volumes, problem dimensions, and query execution modes (Prokopenko et al., 31 Jul 2025).

2. Core search structures and traversal model

The original ArborX design uses a binary BVH in which each node stores an axis-aligned bounding box (AABB) that encloses either two children or exactly one user primitive; in 3D, an AABB is represented by two points, $(x_{\min}, y_{\min}, z_{\min})$ and $(x_{\max}, y_{\max}, z_{\max})$ (Lebrun-Grandié et al., 2019). The 2019 paper describes a fully data-parallel, top-down “linear BVH” build driven by Morton codes, together with iterative stack-based search kernels for range and $k$-nearest-neighbor queries (Lebrun-Grandié et al., 2019). In ray-tracing literature, BVH quality is often discussed through the Surface Area Heuristic,

$$\mathrm{SAH}(\text{node}) = C_{\mathrm{trav}} + \frac{S(A_{\mathrm{left}})}{S(A_{\mathrm{node}})} \cdot C_{\mathrm{left}} + \frac{S(A_{\mathrm{right}})}{S(A_{\mathrm{node}})} \cdot C_{\mathrm{right}},$$

but the original ArborX implementation deliberately chose Morton-code splitting to keep build time low in settings where the tree is rebuilt frequently at each time-step (Lebrun-Grandié et al., 2019).

The exascale paper gives a more implementation-centric account of the same index. It describes a GPU-ready BVH over $N$ primitives, stored in two contiguous arrays of size $2N-1$ in flattened form, with each node holding a $2 \times d$ array of doubles for the AABB, child indices or a rope index after stackless transformation, and leaf payloads as point indices or small lists of point indices (Prokopenko et al., 2024). That work emphasizes rope-based stackless traversal: the right-child pointer of an internal node is replaced with a rope, the index of the next subtree to traverse if the current node fails intersection tests, and every leaf also receives a rope (Prokopenko et al., 2024). Range queries then proceed as purely stackless traversals, while $k$-nearest searches use a per-thread min-heap of size $O(n^2)$0 plus a small stack of nodes still to be visited in distance order (Prokopenko et al., 2024).

ArborX publications use different asymptotic descriptions for BVH construction. The 2019 paper states

$O(n^2)$1

with the parallel sort as the dominant cost (Lebrun-Grandié et al., 2019). The exascale paper describes Morton-code computation in $O(n^2)$2 work, a parallel radix sort on Morton keys in $O(n^2)$3 time with a small constant, Apetrei et al.’s $O(n^2)$4 LBVH construction, and an $O(n^2)$5 rope-recovery pass, summarizing this as $O(n^2)$6 work, $O(n^2)$7 span, memory $O(n^2)$8 (Prokopenko et al., 2024). By contrast, the 2.0 technical report states “BVH build (Apetrei): $O(n^2)$9” and gives per-query traversal as approximately $(x_{\min}, y_{\min}, z_{\min})$0 (Prokopenko et al., 31 Jul 2025). This suggests that the literature is using different accounting conventions for the role of sorting, rather than presenting a single canonical asymptotic statement.

A similar historical variation appears in memory-layout descriptions. The 2019 paper emphasizes Structure-of-Arrays storage for node minima, maxima, and child indices (Lebrun-Grandié et al., 2019). The exascale paper describes an Array-of-Structs-of-Arrays (AoSoA) layout to balance vector-load efficiency with coalescing (Prokopenko et al., 2024). The 2.0 report again describes node data in Structure of Arrays (SoA) for coalesced GPU loads, with leaf pointers to value indices in a separate array (Prokopenko et al., 31 Jul 2025). A plausible implication is that the precise internal layout is an implementation detail that has evolved with hardware tuning and release objectives.

3. Interface evolution and generalized API in version 2.0

The 1.x interface was intentionally narrow. The technical report’s recap of API v1 shows a BVH<MemorySpace> class with a hardwired bounding_volume_type = ArborX::Box, a constructor from primitives, and a query() method returning only indices and offsets in CSR-style form (Prokopenko et al., 31 Jul 2025). Its listed limitations are explicit: hardwired to 3D AABBs, no execution-space argument leading to global fencing, returns only indices plus offsets, and no callbacks or user-data in leaves (Prokopenko et al., 31 Jul 2025).

Version 2.0 replaces that design with a templated BVH<MemorySpace, Value, IndexableGetter, BoundingVolume> and an explicit execution-space argument on both construction and query (Prokopenko et al., 31 Jul 2025). The report identifies the key expansions as arbitrary Value types and a user-supplied IndexableGetter for any geometry, an enum of bounding volumes via template argument, explicit ExecutionSpace for overlap with other Kokkos kernels, and three overloads of query() for zero, one, or many outputs per match (Prokopenko et al., 31 Jul 2025). The class also exposes size(), empty(), and bounds() (Prokopenko et al., 31 Jul 2025).

The architectural rationale is tightly tied to Kokkos. In version 2.0, all data lives in Kokkos::View<> containers parameterized by a MemorySpace, and both construction and query routines take a Kokkos execution space so that work can be dispatched on arbitrary GPU streams or CPU threads (Prokopenko et al., 31 Jul 2025). The same report states that vendor-tuned kernels for sorting and reductions use Thrust, rocThrust, and oneDPL (Prokopenko et al., 31 Jul 2025). The generalized interface is summarized by three design points that materially alter user integration: arbitrary user data types, flexible bounding volumes, and dimension up to 10 on both CPU and GPU (Prokopenko et al., 31 Jul 2025).

This interface evolution also clarifies a recurring misunderstanding about ArborX output semantics. In 1.x, the library returned only materialized search results in CSR-style arrays (Prokopenko et al., 31 Jul 2025). In 2.0, callbacks are a first-class execution model, so result processing can happen in place during traversal and large temporary arrays need not be materialized at all (Prokopenko et al., 31 Jul 2025).

4. Query modes, callbacks, and algorithmic expansion

Callbacks are the main semantic addition in ArborX 2.0. The technical report defines them as a mechanism that lets users run arbitrary code on each (predicate, value) match inside traversal (Prokopenko et al., 31 Jul 2025). It presents two forms. The “pure callback” variant executes user code for each match without stored output, while the “callback with output” variant requires output_type and uses an emitter to append user-defined outputs (Prokopenko et al., 31 Jul 2025). The stated purpose is to allow in-place processing of query results without materializing large temporary arrays (Prokopenko et al., 31 Jul 2025). Traversal can also terminate early when callbacks return a Terminate tag (Prokopenko et al., 31 Jul 2025).

ArborX 2.0 is no longer restricted to a single search data structure. The report adds a brute-force index for cases where “for small $(x_{\min}, y_{\min}, z_{\min})$1 or degenerate point clouds, a brute-force index is often competitive” (Prokopenko et al., 31 Jul 2025). Its internal organization is simple—values in a contiguous Kokkos::View<Value*> and a nested Kokkos::parallel_for over queries and data—and its complexity is given as trivial build $(x_{\min}, y_{\min}, z_{\min})$2 and query

$(x_{\min}, y_{\min}, z_{\min})$3

for $(x_{\min}, y_{\min}, z_{\min})$4 queries and $(x_{\min}, y_{\min}, z_{\min})$5 primitives (Prokopenko et al., 31 Jul 2025). The same document contrasts this with BVH behavior as $(x_{\min}, y_{\min}, z_{\min})$6 build and $(x_{\min}, y_{\min}, z_{\min})$7 query (Prokopenko et al., 31 Jul 2025).

The algorithmic envelope has also widened. The 2.0 release lists ray tracing predicates, density-based clustering (DBSCAN), EMST, and MLS interpolation among supported algorithms (Prokopenko et al., 31 Jul 2025). For ray tracing, geometry is represented as Ray = (origin, direction), and predicates include nearest(k) for the first $(x_{\min}, y_{\min}, z_{\min})$8 intersections, intersect for all intersections, and ordered_intersect for hits sorted by distance (Prokopenko et al., 31 Jul 2025). The report describes a traversal pipeline of ray–node AABB tests, leaf-level ray–geometry intersection for triangle, box, or sphere, and early termination or sorting buffers depending on predicate semantics (Prokopenko et al., 31 Jul 2025).

For DBSCAN, the report gives the standard formulation with points $(x_{\min}, y_{\min}, z_{\min})$9, parameters $(x_{\max}, y_{\max}, z_{\max})$0 and minPts, neighborhood

$(x_{\max}, y_{\max}, z_{\max})$1

and core points determined by $(x_{\max}, y_{\max}, z_{\max})$2 (Prokopenko et al., 31 Jul 2025). It lists two implementations—FDBSCAN for sparse data and FDBSCAN-DenseBox for high local density—and states that both use ArborX for $(x_{\max}, y_{\max}, z_{\max})$3-neighborhood queries in parallel (Prokopenko et al., 31 Jul 2025). EMST is described as a complete graph with weights $(x_{\max}, y_{\max}, z_{\max})$4 and Kruskal’s or Prim’s algorithm accelerated by nearest-neighbor search, while MLS interpolation uses nearest queries to define local windows for weighted least squares (Prokopenko et al., 31 Jul 2025).

5. Distributed search and exascale DBSCAN workflows

Distributed search in ArborX 2.0 is organized through ArborX::DistributedTree, which wraps a local BVH together with a global BVH over local AABBs (Prokopenko et al., 31 Jul 2025). Each MPI rank builds a local $(x_{\max}, y_{\max}, z_{\max})$5 over its data, all ranks exchange local AABBs to build a coarse global $(x_{\max}, y_{\max}, z_{\max})$6, and each query is first routed through $(x_{\max}, y_{\max}, z_{\max})$7 to a subset of ranks before local refinement on the destination ranks (Prokopenko et al., 31 Jul 2025). The report gives per-query communication complexity as $(x_{\max}, y_{\max}, z_{\max})$8 messages plus volume proportional to the number of overlapping ranks, and local cost as $(x_{\max}, y_{\max}, z_{\max})$9 (Prokopenko et al., 31 Jul 2025). GPU-aware MPI support is explicitly part of the 2.0 design goals (Prokopenko et al., 31 Jul 2025).

The exascale paper provides the most detailed distributed DBSCAN formulation. Within each MPI rank, all $k$0 points are inserted into a local BVH, and boundary points are communicated in a halo of radius $k$1 to neighbors so that global connectivity is preserved (Prokopenko et al., 2024). The paper states that this ghost exchange costs $k$2 messages in a regular domain decomposition and later summarizes MPI exchange as $k$3 data volume in one neighbor-round (Prokopenko et al., 2024). On the GPU, the computation is split into two Kokkos kernels: core identification, which issues a range query with early termination once the count hits minPts, and cluster merging, which performs pairwise traversal over all $k$4 with $k$5 and atomically unions points in a disjoint-set structure whenever thread $k$6 visits $k$7 and $k$8 is core (Prokopenko et al., 2024). Because only pairs with $k$9 are processed, every edge is handled exactly once (Prokopenko et al., 2024).

The complexity statement for this exascale DBSCAN is explicit. Build BVH is $$\mathrm{SAH}(\text{node}) = C_{\mathrm{trav}} + \frac{S(A_{\mathrm{left}})}{S(A_{\mathrm{node}})} \cdot C_{\mathrm{left}} + \frac{S(A_{\mathrm{right}})}{S(A_{\mathrm{node}})} \cdot C_{\mathrm{right}},$$0; the core phase is $$\mathrm{SAH}(\text{node}) = C_{\mathrm{trav}} + \frac{S(A_{\mathrm{left}})}{S(A_{\mathrm{node}})} \cdot C_{\mathrm{left}} + \frac{S(A_{\mathrm{right}})}{S(A_{\mathrm{node}})} \cdot C_{\mathrm{right}},$$1 worst case but practically $$\mathrm{SAH}(\text{node}) = C_{\mathrm{trav}} + \frac{S(A_{\mathrm{left}})}{S(A_{\mathrm{node}})} \cdot C_{\mathrm{left}} + \frac{S(A_{\mathrm{right}})}{S(A_{\mathrm{node}})} \cdot C_{\mathrm{right}},$$2 with early stops; the merge phase is $$\mathrm{SAH}(\text{node}) = C_{\mathrm{trav}} + \frac{S(A_{\mathrm{left}})}{S(A_{\mathrm{node}})} \cdot C_{\mathrm{left}} + \frac{S(A_{\mathrm{right}})}{S(A_{\mathrm{node}})} \cdot C_{\mathrm{right}},$$3, where $$\mathrm{SAH}(\text{node}) = C_{\mathrm{trav}} + \frac{S(A_{\mathrm{left}})}{S(A_{\mathrm{node}})} \cdot C_{\mathrm{left}} + \frac{S(A_{\mathrm{right}})}{S(A_{\mathrm{node}})} \cdot C_{\mathrm{right}},$$4 is the number of edges within $$\mathrm{SAH}(\text{node}) = C_{\mathrm{trav}} + \frac{S(A_{\mathrm{left}})}{S(A_{\mathrm{node}})} \cdot C_{\mathrm{left}} + \frac{S(A_{\mathrm{right}})}{S(A_{\mathrm{node}})} \cdot C_{\mathrm{right}},$$5 in rank $$\mathrm{SAH}(\text{node}) = C_{\mathrm{trav}} + \frac{S(A_{\mathrm{left}})}{S(A_{\mathrm{node}})} \cdot C_{\mathrm{left}} + \frac{S(A_{\mathrm{right}})}{S(A_{\mathrm{node}})} \cdot C_{\mathrm{right}},$$6 and $$\mathrm{SAH}(\text{node}) = C_{\mathrm{trav}} + \frac{S(A_{\mathrm{left}})}{S(A_{\mathrm{node}})} \cdot C_{\mathrm{left}} + \frac{S(A_{\mathrm{right}})}{S(A_{\mathrm{node}})} \cdot C_{\mathrm{right}},$$7 is the inverse Ackermann function for Union-Find (Prokopenko et al., 2024). On cosmology data at late times, the paper reports $$\mathrm{SAH}(\text{node}) = C_{\mathrm{trav}} + \frac{S(A_{\mathrm{left}})}{S(A_{\mathrm{node}})} \cdot C_{\mathrm{left}} + \frac{S(A_{\mathrm{right}})}{S(A_{\mathrm{node}})} \cdot C_{\mathrm{right}},$$8 (Prokopenko et al., 2024).

The same work ties these methods to concrete cosmology settings. Its benchmark setup uses Summit, with $$\mathrm{SAH}(\text{node}) = C_{\mathrm{trav}} + \frac{S(A_{\mathrm{left}})}{S(A_{\mathrm{node}})} \cdot C_{\mathrm{left}} + \frac{S(A_{\mathrm{right}})}{S(A_{\mathrm{node}})} \cdot C_{\mathrm{right}},$$9 nodes and $N$0 V100 GPUs per node, and a HACC gravity + hydro box of side $N$1 with $N$2 particles, sampled per rank to $N$3 dark-matter points (Prokopenko et al., 2024). The clustering radius is parameterized by

$N$4

giving $N$5 (Prokopenko et al., 2024). The paper further reports that HACC runs $N$6 gravity-only on $N$7 GPUs for $N$8 long-range steps, invoking DBSCAN on each step, and that hydrodynamic runs with $N$9 particles can perform substructure finding every step in situ (Prokopenko et al., 2024).

A major practical consequence is the migration of halo and galaxy finding from offline post-processing to in-situ analysis. The exascale paper states that, without ArborX, HACC reserved DBSCAN for offline post-processing, leading to $2N-1$0 of dumped particle data (Prokopenko et al., 2024). With ArborX, cluster catalogs and galaxy identifications, with minPts=10 for stars, are produced on the fly at each major timestep with negligible disruption to the solver (Prokopenko et al., 2024).

6. Performance characteristics and empirical evaluation

The 2019 evaluation establishes ArborX as competitive with CPU-only state-of-the-art libraries even before the exascale additions. In single-threaded comparisons on the Elseberg synthetic data sets, the paper reports ArborX build at approximately $2N-1$1 faster than nanoflann and on par with Boost.Geometry.Index for $2N-1$2 (Lebrun-Grandié et al., 2019). For $2N-1$3NN query throughput on filled cases, it reports ArborX as $2N-1$4–$2N-1$5 faster for $2N-1$6, with even higher gains on hollow cases (Lebrun-Grandié et al., 2019). For range queries, the paper gives $2N-1$7–$2N-1$8 higher throughput in the 2-pass mode, about another $2N-1$9 improvement for the 1-pass mode when the buffer estimate holds, and up to $2 \times d$0–$2 \times d$1 higher throughput for hollow cases dominated by empty queries (Lebrun-Grandié et al., 2019). In terms of spatial search rate, it reports roughly $2 \times d$2 million boxes per second for filled cases and roughly $2 \times d$3 million boxes per second for hollow cases (Lebrun-Grandié et al., 2019). On Summit, one V100 GPU is reported as $2 \times d$4–$2 \times d$5 faster than $2 \times d$6 POWER9 cores at SMT4 for $2 \times d$7, with range and $2 \times d$8NN queries achieving more than $2 \times d$9 billion node-tests per second on V100 (Lebrun-Grandié et al., 2019).

The exascale paper isolates a timeline of DBSCAN-specific performance improvements across ArborX releases. It attributes a $k$0 time reduction and halved memory overhead to FDbscan with callbacks and early termination in v1.4; a $k$1 traversal-cost reduction to Apetrei build plus rope-stackless traversal in v1.6–v1.7; a further $k$2 speedup on clustered data to 64-bit Morton codes in v1.9; and another $k$3 gain to pairwise stackless traversal in v1.10 (Prokopenko et al., 2024). Cumulatively, the paper reports per-rank DBSCAN runtime dropping from approximately $k$4 to under $k$5, a $k$6 speedup (Prokopenko et al., 2024). For full HACC workflows, it reports ArborX-GPU DBSCAN at $k$7–$k$8 speedup over OpenMP-CPU and an approximately $k$9 end-to-end acceleration of the production code (Prokopenko et al., 2024).

The 2.0 technical report complements those release-history numbers with representative build-and-query trends across architectures. For strong scaling on an NVIDIA A100 with $O(n^2)$00 primitives and $O(n^2)$01 queries, it reports build time improving from $O(n^2)$02 on a $O(n^2)$03-core AVX2 CPU to $O(n^2)$04 on A100, a $O(n^2)$05 speedup, and intersect-query time improving from $O(n^2)$06 to $O(n^2)$07, a $O(n^2)$08 speedup (Prokopenko et al., 31 Jul 2025). For weak scaling in distributed mode with fixed per-rank $O(n^2)$09 and $O(n^2)$10 over $O(n^2)$11 MPI ranks, the report states that total time is approximately constant build plus $O(n^2)$12 query routing (Prokopenko et al., 31 Jul 2025). Its GPU breakdown for $O(n^2)$13 boxes and $O(n^2)$14 points reports $O(n^2)$15 for Morton sort, $O(n^2)$16 for BVH build, $O(n^2)$17 for spatial query, and $O(n^2)$18 for ray tracing with $O(n^2)$19 (Prokopenko et al., 31 Jul 2025). The same report describes a log–log plot of query time versus $O(n^2)$20 whose slopes match the theoretical exponents for $O(n^2)$21 build and near-$O(n^2)$22 query (Prokopenko et al., 31 Jul 2025).

Taken together, these evaluations define ArborX as a portable HPC search substrate whose performance story is not limited to raw BVH traversal. The reported gains depend on a sequence of interlocking mechanisms: Kokkos-based back-end portability, 64-bit Morton codes, Apetrei-style build strategies, stackless or pairwise traversals, early-out callbacks, distributed routing through a two-level index, and algorithm-specific integrations such as DBSCAN and ray tracing (Prokopenko et al., 2024, Prokopenko et al., 31 Jul 2025).