P³T Scheme: Hybrid N-body Simulation Method

• P³T scheme is a hybrid method that merges O(N²) particle-particle and O(N log N) particle-tree approaches to simulate collisional gravitational systems.
• It employs Hamiltonian splitting and robust free-fall-based switching criteria to ensure sharp energy conservation across different dynamical regimes.
• By integrating Hermite time-stepping with adaptive multipole tree dynamics, the scheme offers a versatile balance of accuracy and computational efficiency.

The P$^3$T scheme (particle–particle–particle–tree) is a hybrid computational method designed for high-precision simulations of collisional gravitational $N$-body systems, such as star clusters. It achieves both scalability and accuracy by combining a $O(N\log N)$ particle-tree (PT) algorithm for long-range forces with a direct $O(N^2)$ particle-particle (PP) solver for short-range interactions. Key contributions include the introduction of robust switching criteria between the PT and PP regions, especially a free-fall-based criterion, that enables efficient and accurate energy conservation across regimes with varying velocity dispersion and initial conditions (Wang et al., 12 Jan 2026).

1. Hamiltonian Splitting and Hybrid Dynamics

In the P$^3$T formalism, the total Hamiltonian of an $N$-body system,

$$H = \sum_{i=1}^N \frac{p_i^2}{2m_i} - \sum_{i<j}^{N} \frac{G m_i m_j}{r_{ij}},$$

is decomposed into two components using a smooth changeover function $W(r_{ij})$:

• The hard part,

$$H_{\mathrm{hard}} = \sum_i \frac{p_i^2}{2m_i} - \sum_{i<j} \frac{G m_i m_j}{r_{ij}} W(r_{ij}),$$

is integrated with a direct Hermite $N^2$ solver, efficiently handling close encounters.

• The soft part,

$N$0

is evaluated by a tree code using a multipole expansion up to quadrupole.

For $N$1, $N$2 so that the direct solver dominates; for $N$3, $N$4 and the tree dominates; in the transition region, contributions blend smoothly (Wang et al., 12 Jan 2026).

2. Tree Algorithm Structure and Force Computation

The particle-tree component is based on a hierarchical Barnes–Hut octree:

• Each node contains multipole moments (up to quadrupole or higher).
• The node-opening condition,

$N$5

controls the acceptance of a node (size $N$6 at center $N$7 relative to its distance $N$8 from the target particle $N$9), with $O(N\log N)$0 an accuracy parameter.

• Forces are computed at $O(N\log N)$1 scaling per tree step.
• The tree integration uses a shared timestep $O(N\log N)$2, which is set to resolve the shortest relevant two-body free-fall time.

The coupling to the Hermite integrator is achieved by matching the maximum Hermite time block to $O(N\log N)$3 and by using an eighth-order changeover function $O(N\log N)$4 for smooth force derivatives at interface boundaries (Wang et al., 12 Jan 2026).

3. Direct Particle-Particle Solver and Regularization

For close pairwise interactions, the PP region is integrated using a fourth-order Hermite scheme with block timesteps:

• The individual timestep is determined by Aarseth’s formula,

$O(N\log N)$5

where $O(N\log N)$6 is the $O(N\log N)$7-th derivative of acceleration, $O(N\log N)$8, and $O(N\log N)$9.

For compact subsystems (e.g., binaries, higher-order multiples), slow-down algorithmic regularization (SDAR) is triggered. SDAR is based on a time-transformed leapfrog integrator, allowing efficient regularized evolution of tight subsystems while preserving the symplectic structure (Wang et al., 12 Jan 2026).

4. PT-PP Switching Criteria

Effective hybridization requires a robust method to determine which pairs are handled by the tree or by the direct $O(N^2)$0 solver. Two principal switching criteria have been employed:

• Free-fall-based criterion: The soft timestep $O(N^2)$1 is set to resolve a two-body circular orbital period at radius $O(N^2)$2 by at least $O(N^2)$3 steps,

$O(N^2)$4

and

$O(N^2)$5

Typically, $O(N^2)$6 and the maximum relevant $O(N^2)$7 are used to set a global timestep (Wang et al., 12 Jan 2026).

• Velocity-dispersion-based ($O(N^2)$8-based) criterion: Adopted in prior PETAR versions, the switch is defined by

$O(N^2)$9

where $^3$0 is the global 3D velocity dispersion and $^3$1 is a tuning parameter (Wang et al., 12 Jan 2026).

The PP region is defined by all pairs with $^3$2, where $^3$3.

5. Quantitative Performance: Accuracy and Application Regimes

Comparative studies in virial-equilibrium, subvirial, and fractal initial conditions reveal differing regimes of superiority:

• In virial-equilibrium clusters ($^3$4, $^3$5, $^3$6), both criteria predict the same critical $^3$7. For $^3$8, free-fall-based switching yields larger $^3$9 and smaller energy errors (e.g., at $N$0, $N$1 for free-fall vs $N$2 for $N$3-based) (Wang et al., 12 Jan 2026).
• For $N$4, the transition timestep scales as $N$5; thus, $N$6-based switching outperforms in tight, high-$N$7 systems.
• In subvirial or fractal initial conditions, $N$8-based switching fails to maintain energy conservation due to a poor estimate of the dynamical timescale; free-fall-based switching, in contrast, demonstrates monotonic improvement of energy conservation with reduced $N$9 (see Fig. 20, 23 in (Wang et al., 12 Jan 2026)).
• For low-$$H = \sum_{i=1}^N \frac{p_i^2}{2m_i} - \sum_{i<j}^{N} \frac{G m_i m_j}{r_{ij}},$$0 associations ($$H = \sum_{i=1}^N \frac{p_i^2}{2m_i} - \sum_{i<j}^{N} \frac{G m_i m_j}{r_{ij}},$$1) and embedded clusters, the free-fall-based method is preferred; for dense nuclear clusters ($$H = \sum_{i=1}^N \frac{p_i^2}{2m_i} - \sum_{i<j}^{N} \frac{G m_i m_j}{r_{ij}},$$2) and $$H = \sum_{i=1}^N \frac{p_i^2}{2m_i} - \sum_{i<j}^{N} \frac{G m_i m_j}{r_{ij}},$$3, $$H = \sum_{i=1}^N \frac{p_i^2}{2m_i} - \sum_{i<j}^{N} \frac{G m_i m_j}{r_{ij}},$$4-based becomes more efficient at fixed accuracy.

6. Practical Recommendations and Algorithmic Settings

Parameter selection is system dependent:

• For $$H = \sum_{i=1}^N \frac{p_i^2}{2m_i} - \sum_{i<j}^{N} \frac{G m_i m_j}{r_{ij}},$$5 and $$H = \sum_{i=1}^N \frac{p_i^2}{2m_i} - \sum_{i<j}^{N} \frac{G m_i m_j}{r_{ij}},$$6, implement the free-fall criterion with $$H = \sum_{i=1}^N \frac{p_i^2}{2m_i} - \sum_{i<j}^{N} \frac{G m_i m_j}{r_{ij}},$$7 to choose $$H = \sum_{i=1}^N \frac{p_i^2}{2m_i} - \sum_{i<j}^{N} \frac{G m_i m_j}{r_{ij}},$$8:

$$H = \sum_{i=1}^N \frac{p_i^2}{2m_i} - \sum_{i<j}^{N} \frac{G m_i m_j}{r_{ij}},$$9

with $W(r_{ij})$0.

• For high-$W(r_{ij})$1, high-$W(r_{ij})$2 clusters (e.g., globular, nuclear star clusters), optionally switch to the $W(r_{ij})$3-based criterion with $W(r_{ij})$4.
• In mixed or evolving environments, the free-fall criterion remains robust as it does not rely on a global $W(r_{ij})$5.
• Barnes-Hut opening angle $W(r_{ij})$6 in $W(r_{ij})$7 balances accuracy and performance.
• For regularization, the PP subsystem threshold $W(r_{ij})$8 ensures efficient treatment of hard binaries (Wang et al., 12 Jan 2026).

A hybrid approach choosing the maximum $W(r_{ij})$9 from both criteria can provide seamless performance across parameter regimes.

7. Advantages, Limitations, and Scope of Application

The P$$H_{\mathrm{hard}} = \sum_i \frac{p_i^2}{2m_i} - \sum_{i<j} \frac{G m_i m_j}{r_{ij}} W(r_{ij}),$$0T scheme, particularly with a free-fall-based switching criterion, features:

• Avoidance of global $$H_{\mathrm{hard}} = \sum_i \frac{p_i^2}{2m_i} - \sum_{i<j} \frac{G m_i m_j}{r_{ij}} W(r_{ij}),$$1 estimation, reducing bias in low-$$H_{\mathrm{hard}} = \sum_i \frac{p_i^2}{2m_i} - \sum_{i<j} \frac{G m_i m_j}{r_{ij}} W(r_{ij}),$$2 and clumpy systems.
• Adaptive resilience, guaranteeing at least $$H_{\mathrm{hard}} = \sum_i \frac{p_i^2}{2m_i} - \sum_{i<j} \frac{G m_i m_j}{r_{ij}} W(r_{ij}),$$3 leapfrog steps per minimal orbit at $$H_{\mathrm{hard}} = \sum_i \frac{p_i^2}{2m_i} - \sum_{i<j} \frac{G m_i m_j}{r_{ij}} W(r_{ij}),$$4 and systematically improving energy conservation as the timestep decreases.
• Suitability for simulations of open/embedded clusters, early-stage star formation regions, and systems with evolving velocity structure (Wang et al., 12 Jan 2026).

Limitations include less efficiency for high-$$H_{\mathrm{hard}} = \sum_i \frac{p_i^2}{2m_i} - \sum_{i<j} \frac{G m_i m_j}{r_{ij}} W(r_{ij}),$$5, high-$$H_{\mathrm{hard}} = \sum_i \frac{p_i^2}{2m_i} - \sum_{i<j} \frac{G m_i m_j}{r_{ij}} W(r_{ij}),$$6 systems, where the $$H_{\mathrm{hard}} = \sum_i \frac{p_i^2}{2m_i} - \sum_{i<j} \frac{G m_i m_j}{r_{ij}} W(r_{ij}),$$7-based criterion may deliver superior computational cost at equivalent accuracy, and the need for empirical correction factors near $$H_{\mathrm{hard}} = \sum_i \frac{p_i^2}{2m_i} - \sum_{i<j} \frac{G m_i m_j}{r_{ij}} W(r_{ij}),$$8 orbits. In applications requiring both scalable long-range force computation and high-accuracy integration of close encounters, the P$$H_{\mathrm{hard}} = \sum_i \frac{p_i^2}{2m_i} - \sum_{i<j} \frac{G m_i m_j}{r_{ij}} W(r_{ij}),$$9T formalism with combined switching criteria offers robust, versatile modeling capabilities.