Phantom N-Body Framework Overview

• Phantom N-Body Framework is a collection of diverse computational methods unified by hardware-aware auto-tuning and minimal manual intervention across various physical domains.
• It integrates hybrid treecode–FMM, distributed Barnes–Hut, SIMD-accelerated kernels, and AMR-based solvers to efficiently handle long-range interactions and complex dynamics.
• Recent developments extend the framework for embedded stellar clusters with advanced integrators and for quantum systems by revealing phantom vortices in coherence space.

The designation Phantom N-Body Framework is applied in the literature to several distinct, technically unrelated research programs rather than to a single canonical software stack. In computational astrophysics and high-performance computing, it denotes black-box hybrid treecode–FMM designs for heterogeneous systems, distributed Barnes–Hut realizations on MPI clusters, SIMD-accelerated collisionless force libraries, and AMR-based Milgromian gravity solvers; in recent star-cluster work it also names a hybrid SPH-plus-collisional-dynamics framework built inside Phantom; and, in a separate quantum many-body context, “phantom” denotes hidden coherence-space vorticity in rotating bosonic systems (Yokota et al., 2011, Brandt, 2022, Tanikawa et al., 2012, Lüghausen et al., 2014, Nagesh et al., 2021, Bernard et al., 7 Aug 2025, Weiner et al., 2014).

1. Scope of the designation

The literature associates the phrase with several frameworks that share neither a single numerical method nor a single physical domain. This suggests that Phantom N-Body Framework functions as a family resemblance term rather than a uniquely defined standard.

Framework line	Core setting	Characteristic idea
Hybrid treecode–FMM	Heterogeneous HPC	Auto-tuned switching among P2P, cell-particle, and cell-cell
Distributed Barnes–Hut	MPI clusters	Hashed octree, branch-node exchange, remote child fetch
Phantom-GRAPE	x86 SIMD	AVX-accelerated pairwise kernels for collisionless solvers
Phantom of RAMSES	AMR galaxy simulations	QUMOND solver and phantom density on RAMSES meshes
DAWN / Phantom N-Body	Embedded clusters	FSI + SDAR inside Phantom with gas, sinks, and H II feedback
Phantom vortices	Rotating bosonic many-body systems	Quantized vorticity in coherence space rather than density

A common misconception is that the term refers to one software package. The published record instead shows multiple uses tied to different scientific problems. Another misconception is that “phantom” always denotes effective dark matter or hidden degrees of freedom. In the MOND context it refers to the effective source term $\rho_{\rm ph}$; in the bosonic many-body context it refers to vortices detectable in coherence but not in the density; and in the DAWN framework it is inherited from the host hydrodynamics code Phantom rather than from a specific physical approximation (Lüghausen et al., 2014, Weiner et al., 2014, Bernard et al., 7 Aug 2025).

2. Black-box fast N-body algorithms on heterogeneous systems

One major usage is the black-box, auto-tuned hybrid treecode–FMM framework for heterogeneous hardware. Its target class is long-range interactions governed by Laplace-type kernels, with stated extensibility to Helmholtz and Stokes potentials via appropriate multipole/local operators. “Black-box” means that the framework automatically chooses algorithmic paths, operator variants, and device scheduling to meet accuracy and performance targets on CPUs, GPUs, and accelerators without requiring manual tuning of internal parameters such as the order $p$, opening angle $\theta$, particles per leaf $n_{\rm leaf}$, or the operator mix (Yokota et al., 2011).

The mathematical core combines hierarchical spatial decomposition with multipole and local expansions. The acceptance condition is expressed through the multipole acceptance criterion

$\theta = \frac{r_t + r_s}{R},$

with interaction accepted if $\theta \le \theta_{\max}$. Within a single dual-tree traversal, the framework may evaluate a cell pair by direct particle–particle interaction, by cell-particle interaction in a treecode style, or by cell-cell interaction in an FMM style. The cost model is written as

$$C_{\rm total} = \alpha N_{\rm P2P} + \beta N_{\rm M2L} + \gamma N_{\rm M2M} + \delta N_{\rm P2M} + \epsilon N_{\rm L2L} + \zeta N_{\rm L2P},$$

where the coefficients are device-specific per-invocation costs measured or updated online. This cost-based switching is the core mechanism by which the same code path can favor cell-cell interactions on CPUs and shift work toward P2P on GPUs when arithmetic intensity and bandwidth behavior make that favorable (Yokota et al., 2011).

The implementation uses adaptive octrees or quadtrees, linear octrees keyed by Morton order, structure-of-arrays particle layouts, compact coefficient buffers, and asynchronous CPU–GPU scheduling. The stated complexity remains the standard one: treecodes scale as $O(N\log N)$ and FMM as $O(N)$. Reported performance figures include a Laplace potential-plus-force run for $N=10^7$ with $p$0 completing in approximately 23 seconds on one NVIDIA GPU in single precision, and weak scaling with 72% efficiency up to 2048 GPUs for an FMM-based vortex method. In this line of work, the phrase Phantom-like chiefly denotes portability, hardware-aware auto-tuning, and minimal user intervention rather than a specific gravitational or astrophysical model (Yokota et al., 2011).

3. Distributed Barnes–Hut realizations

A second line of work treats a scalable distributed Barnes–Hut gravitational framework. The starting point is the standard complexity contrast: direct particle–particle summation scales as $p$1 per step, whereas Barnes–Hut reduces this to $p$2 by aggregating distant particles through a spatial tree. The paper distinguishes two complementary MPI designs: a costzones-based global-tree approach, which builds local trees and merges them into a global tree, and a hashed-octree approach, which keeps trees local, shares only branch nodes, requests remote child cells on demand, and overlaps communication with computation (Brandt, 2022).

The Barnes–Hut force law is given in both unsoftened and softened form, with softening introduced to avoid singularities and stabilize close encounters. Traversal follows the standard opening criterion $p$3, with $p$4 described as a robust default. The framework supports monopole approximations and outlines quadrupole coefficients up to order 2. Time integration uses leapfrog in velocity-Verlet KDK form with fixed $p$5 in the reported experiments (Brandt, 2022).

The major algorithmic distinction is communicational. The global-tree variant relies on gathers, scatters, reductions, and broadcasts, though it is improved by distributed sort and neighbor rebalancing. The hashed-octree variant avoids global tree broadcast: ranks exchange lightweight branch-node sets, then during force evaluation issue nonblocking child requests keyed by Morton or octal identifiers, maintaining a walk list and a defer list to hide latency. The paper reports that the hashed-octree variant, identified as Algorithm 7, delivers the best scalability and even superlinear speedups at low-to-moderate process counts because it avoids global tree construction. Across Barnes–Hut variants, the total energy change after 500 steps is reported as approximately 0.13–0.15%, and parallel results match serial results in that metric. The broader implication is that, for multi-node MPI clusters, the decisive optimization is often communication avoidance rather than further refinement of the basic Barnes–Hut approximation itself (Brandt, 2022).

4. SIMD collisionless kernels and Phantom-GRAPE

A third usage is Phantom-GRAPE, a numerical software library for collisionless N-body simulations on x86 CPUs using SIMD vector instructions, specifically AVX as an enhancement over SSE. Its role is narrower than that of a full tree or particle–mesh code: it accelerates the pairwise interaction kernel for Newtonian gravity and for arbitrary central forces with finite cutoff radius $p$6, making it especially suitable for tree and TreePM solvers with individual timesteps and irregular interaction lists (Tanikawa et al., 2012).

The implementation is built around 256-bit YMM registers, single-precision arithmetic, careful data alignment, and a batching strategy that computes forces from two $p$7-particles on four $p$8-particles simultaneously. For softened Newtonian gravity it evaluates

$p$9

while the arbitrary-force mode tabulates $\theta$0 and interpolates from a cache-resident lookup table. The library uses the approximate reciprocal square root instruction vrsqrtps, with statistical debiasing at initialization and no Newton–Raphson refinement, explicitly prioritizing throughput for collisionless applications (Tanikawa et al., 2012).

The reported performance on an Intel Core i7-2600 processor is concrete: for Newtonian forces, $\theta$1 interactions per second on one core and $\theta$2 on four cores; for arbitrary central forces, $\theta$3 on one core and $\theta$4 on four cores. The paper emphasizes that performance depends weakly on the number of particles, unlike typical GPU behavior. Relative force errors are reported as typically $\theta$5 for Newtonian gravity and $\theta$6 for the S2-based PPPM/TreePM kernel. This library therefore occupies a specific niche: it is not a complete “phantom framework” in the architectural sense of distributed or AMR codes, but a low-overhead SIMD kernel layer that complements hierarchical solvers in CPU-centric and MPI-dominated environments (Tanikawa et al., 2012).

5. AMR Milgromian gravity and Phantom of RAMSES

In galactic dynamics, the best-known phantom terminology is Phantom of RAMSES (PoR), a customized version of RAMSES that implements QUMOND on an adaptive mesh refinement hierarchy while retaining RAMSES’ particle, gas, and hydrodynamic infrastructure. The defining equations are the Newtonian Poisson solve

$\theta$7

followed by the QUMOND solve

$\theta$8

or, equivalently,

$\theta$9

Here $n_{\rm leaf}$0 is the effective phantom dark matter density, a Newtonian reinterpretation of the MOND enhancement sourced entirely by baryons (Lüghausen et al., 2014).

PoR uses RAMSES multigrid Gauss–Seidel relaxation, isolated boundary conditions, and a “simple” QUMOND interpolation function. The later user guide specifies that PoR is a patch to the 2015 version of RAMSES and documents the practical run-time controls: mond=.true. or .false. to toggle MOND and Newtonian gravity, a0_ms2=1.2e-10, Activate_g_ext for the External Field Effect option in the hydro patch, and patch variants such as patch/phantom_staticparts and patch/hydro/phantom_extfield. It also documents DICE-based initial-condition generation for isolated and interacting disc galaxies and post-processing tools such as extract_por and extract_por_sfr (Nagesh et al., 2021).

The scientific demonstrations are equally central to the framework’s identity. Static tests recover the expected deep-MOND acceleration and phantom-density profiles for point-mass and Plummer configurations. Dynamical tests show stable spherical equilibria over more than 100 crossing times. Disk simulations produce rapid outer spiral substructure, strong bar formation by approximately 1 Gyr, nearly constant pattern speed over several Gyr, and a phantom density distribution comprising a near-isothermal oblate halo together with a phantom dark disk aligned with the stellar disk. A recurring contrast with Newtonian-plus-dark-matter runs is that PoR allows side-by-side comparisons on identical AMR meshes and with identical baryonic initial conditions, isolating the consequences of the gravity law itself (Lüghausen et al., 2014).

6. Embedded stellar clusters in Phantom: FSI, SDAR, and feedback

A more recent and substantially different use of the designation appears in the DAWN programme, which introduces a new Phantom N-Body framework for the formation and early evolution of stellar clusters inside a molecular cloud environment. The framework is built directly into the SPH code Phantom and is motivated by the coexistence of collisional stellar dynamics, turbulent gas collapse, and feedback over a large dynamical range. The original difficulty was that standard SPH runs in Phantom using second-order leapfrog struggled to integrate long, collisional cluster evolutions without prohibitive timesteps and accumulating truncation error (Bernard et al., 7 Aug 2025).

The algorithmic solution is twofold. First, Phantom’s slow–fast RESPA force split is retained, but the fast-level leapfrog is replaced by a fourth-order forward symplectic integrator (FSI) following Dehnen & Hernandez and Rantala et al. Its modified mid-kick uses the gradient acceleration and Omelyan extrapolation. Second, close encounters and hard binaries are handled through algorithmic regularization coupled to slowdown (SDAR). Subgroups are detected at each $n_{\rm leaf}$1 using neighbors within $n_{\rm leaf}$2 pc plus an imminent-fly-by criterion with $n_{\rm leaf}$3. These subgroups are then integrated with a time-transformed leapfrog and a 6th- or 8th-order Yoshida symplectic scheme, while hard binaries can be slowed through

$n_{\rm leaf}$4

with adaptive $n_{\rm leaf}$5 modulated by the Mikkola–Aarseth tidal measure and limited by $n_{\rm leaf}$6 when perturbations vary too rapidly (Bernard et al., 7 Aug 2025).

The surrounding physical model is equally specific. Star–star, sink–sink, and sink–star forces are unsoftened; star–gas and sink–gas interactions are softened, with $n_{\rm leaf}$7 au for stars in star–gas interactions. Star formation proceeds through sinks created at $n_{\rm leaf}$8 with $n_{\rm leaf}$9 au, subject to the binding check $\theta = \frac{r_t + r_s}{R},$0. Each sink accretes for $\theta = \frac{r_t + r_s}{R},$1 Myr, forms 1–5 star seeds after $\theta = \frac{r_t + r_s}{R},$2 Myr, and then dissolves into a bound, virialized mini-cluster whose stellar masses satisfy $\theta = \frac{r_t + r_s}{R},$3. Massive stars above $\theta = \frac{r_t + r_s}{R},$4 drive H II region expansion through a particle-tagging photon-budget scheme with

$\theta = \frac{r_t + r_s}{R},$5

The fiducial GMC initial conditions are $\theta = \frac{r_t + r_s}{R},$6, $\theta = \frac{r_t + r_s}{R},$7 pc, $\theta = \frac{r_t + r_s}{R},$8 K, $\theta = \frac{r_t + r_s}{R},$9, $\theta \le \theta_{\max}$0, and 3.5 million SPH particles (Bernard et al., 7 Aug 2025).

The reported testcase simulations ECF1, ECF2, and ECF3 show that the placement of the first massive stars and the coherence of multiple H II regions determine whether the cloud undergoes rapid dismantling or a second, low-mass-dominated collapse. Cloud-scale star formation efficiencies are 14–18%, the mean star formation rates are $\theta \le \theta_{\max}$1 for ECF1 and ECF2 and $\theta \le \theta_{\max}$2 for ECF3, and the final star count in ECF2 reaches $\theta \le \theta_{\max}$3. On the numerical side, the paper reports that restarting near the birth of a hard binary without SDAR required 48 hours to complete only one-fifth of the usual step-to-dump interval, whereas with FSI+SDAR the full interval took 3.1 hours, a speed-up of approximately 80. The main late-time bottleneck is no longer the N-body integrator itself but the direct $\theta \le \theta_{\max}$4 star/sink–gas force calculation, which accounts for more than 80% of walltime once $\theta \le \theta_{\max}$5 (Bernard et al., 7 Aug 2025).

7. Quantum many-body usage: phantom vortices in coherence space

The broadest extension of the terminology appears in rotating interacting bosonic systems, where the “phantom” label no longer denotes a gravitational framework at all but a many-body diagnostic of hidden angular momentum. In that setting, a phantom vortex is defined as a quantized phase defect that resides in the many-body coherence of a single-component Bose gas, not in its particle density. It is visible as a topological defect in the first-order coherence function $\theta \le \theta_{\max}$6 and in the phases of natural orbitals extracted from the one-body reduced density matrix, but it does not appear as a conventional density core. Angular momentum is therefore hidden from density-based diagnostics (Weiner et al., 2014).

The formalism uses MCTDHB, in which the many-body state is expanded over $\theta \le \theta_{\max}$7 time-dependent variational orbitals and the Gross–Pitaevskii description is recovered only in the special case $\theta \le \theta_{\max}$8. The one-body density matrix

$\theta \le \theta_{\max}$9

is diagonalized to obtain natural orbitals $$C_{\rm total} = \alpha N_{\rm P2P} + \beta N_{\rm M2L} + \gamma N_{\rm M2M} + \delta N_{\rm P2M} + \epsilon N_{\rm L2L} + \zeta N_{\rm L2P},$$0 and occupations $$C_{\rm total} = \alpha N_{\rm P2P} + \beta N_{\rm M2L} + \gamma N_{\rm M2M} + \delta N_{\rm P2M} + \epsilon N_{\rm L2L} + \zeta N_{\rm L2P},$$1, and the total angular momentum is partitioned as

$$C_{\rm total} = \alpha N_{\rm P2P} + \beta N_{\rm M2L} + \gamma N_{\rm M2M} + \delta N_{\rm P2M} + \epsilon N_{\rm L2L} + \zeta N_{\rm L2P},$$2

Phantom vortices are identified either through zeros and $$C_{\rm total} = \alpha N_{\rm P2P} + \beta N_{\rm M2L} + \gamma N_{\rm M2M} + \delta N_{\rm P2M} + \epsilon N_{\rm L2L} + \zeta N_{\rm L2P},$$3 phase windings in the natural orbitals or through spots where $$C_{\rm total} = \alpha N_{\rm P2P} + \beta N_{\rm M2L} + \gamma N_{\rm M2M} + \delta N_{\rm P2M} + \epsilon N_{\rm L2L} + \zeta N_{\rm L2P},$$4 and the coherence phase winds by $$C_{\rm total} = \alpha N_{\rm P2P} + \beta N_{\rm M2L} + \gamma N_{\rm M2M} + \delta N_{\rm P2M} + \epsilon N_{\rm L2L} + \zeta N_{\rm L2P},$$5 (Weiner et al., 2014).

The representative simulations involve $$C_{\rm total} = \alpha N_{\rm P2P} + \beta N_{\rm M2L} + \gamma N_{\rm M2M} + \delta N_{\rm P2M} + \epsilon N_{\rm L2L} + \zeta N_{\rm L2P},$$6 bosons in a rotating anisotropic 2D trap with $$C_{\rm total} = \alpha N_{\rm P2P} + \beta N_{\rm M2L} + \gamma N_{\rm M2M} + \delta N_{\rm P2M} + \epsilon N_{\rm L2L} + \zeta N_{\rm L2P},$$7 orbitals, weak interactions, and zero temperature. Fragmentation begins around $$C_{\rm total} = \alpha N_{\rm P2P} + \beta N_{\rm M2L} + \gamma N_{\rm M2M} + \delta N_{\rm P2M} + \epsilon N_{\rm L2L} + \zeta N_{\rm L2P},$$8, the occupations approach a completely fragmented configuration by $$C_{\rm total} = \alpha N_{\rm P2P} + \beta N_{\rm M2L} + \gamma N_{\rm M2M} + \delta N_{\rm P2M} + \epsilon N_{\rm L2L} + \zeta N_{\rm L2P},$$9, and at $O(N\log N)$0 a typical occupation distribution is approximately 40.8%, 25.0%, 20.8%, and 13.5% across the four natural orbitals. The many-body angular momentum stabilizes at approximately $O(N\log N)$1 per particle even though the density shows no visible vortices, whereas under identical parameters the GP calculation yields less than $O(N\log N)$2 per particle and no vortices. The paper identifies two many-body nucleation mechanisms—node mutation and slow orbital–orbital vorticity transfer—and argues that coherence-based interferometric measurements of $O(N\log N)$3 would be the natural route to experimental detection. In this context, the “phantom” concept shifts from effective sources or inherited code names to hidden topology in coherence space, but the recurring theme remains the same: physically relevant structure is present even when conventional diagnostics appear featureless (Weiner et al., 2014).