MSTAR AR-CHAIN Integrator

• MSTAR AR-CHAIN integrator is a high-precision, massively parallel algorithm for gravitational N-body problems, achieving machine-precision accuracy in close encounters.
• It utilizes a minimum spanning tree-based coordinate construction to significantly reduce numerical errors, outperforming classical chain methods by an order of magnitude.
• The integrator combines Gragg–Bulirsch–Stoer extrapolation with two-level parallelisation, enabling efficient simulations of galactic nuclei and SMBH binaries with high particle counts.

The MSTAR AR-CHAIN integrator is a high-precision, massively parallel algorithm for integrating gravitational $N$-body systems, specifically engineered to achieve machine-precision accuracy ($|\Delta E / E| \lesssim 10^{-14}$) in configurations dominated by close encounters. Building on the algorithmic regularisation (AR) framework that underlies the AR-CHAIN method, MSTAR introduces a combinatorial coordinate construction based on minimum spanning trees (MST), together with advanced extrapolation schemes and two-level parallelisation that together circumvent computational bottlenecks of earlier approaches. The method enables direct summation integrations of regularised subsystems at a scale ($N_\mathrm{part} \sim 10^4$) and with a level of accuracy previously unattainable in galactic dynamics and stellar cluster modelling, facilitating next-generation simulations of galaxy nuclei and supermassive black hole (SMBH) binaries (Rantala et al., 2020).

1. Algorithmic Regularisation: Hamiltonian Structure and Time Transformation

MSTAR, like AR-CHAIN, operates by restructuring the classical Newtonian $N$-body Hamiltonian

$$H(\{ \mathbf r_i,\mathbf v_i\}) = \sum_{i=1}^N\frac{1}{2}m_i\,\mathbf v_i^2 - \sum_{i<j}\frac{G\,m_i\,m_j}{|\mathbf r_j-\mathbf r_i|},$$

into a regularised form by extending phase space. Physical time $t$ becomes a canonical coordinate with conjugate momentum $B = -H$, and a fictitious regularisation time $s$ is introduced via the transformation

$$\frac{dt}{ds} = \frac{1}{\alpha\,U+\beta\,\Omega+\gamma}.$$

For both MSTAR and AR-CHAIN, the logarithmic Hamiltonian prescription $(\alpha,\beta,\gamma)=(1,0,0)$ yields $|\Delta E / E| \lesssim 10^{-14}$0, eliminating the singular behaviour of $|\Delta E / E| \lesssim 10^{-14}$1 as pairs approach. The transformed Hamiltonian becomes

$|\Delta E / E| \lesssim 10^{-14}$2

permitting a splitting into drift (momentum only) and kick (position only) operators, on which a second-order leapfrog in $|\Delta E / E| \lesssim 10^{-14}$3 is exact for Keplerian orbits. As $|\Delta E / E| \lesssim 10^{-14}$4, $|\Delta E / E| \lesssim 10^{-14}$5 ensures that "physical time stalls at collisions," rendering the equations of motion manifestly regular and free of force divergences (Rantala et al., 2020, Wang et al., 2021).

2. MST-Based Coordinate Construction and Stability Advantages

Unlike AR-CHAIN, which employs a single, sequential chain of coordinate differences between adjacent particles, MSTAR constructs a minimum spanning tree on the particle positions, selecting $|\Delta E / E| \lesssim 10^{-14}$6 edges that connect all particles with minimal total length via (e.g.) Prim's algorithm. Each edge is directed from a child $|\Delta E / E| \lesssim 10^{-14}$7 to a unique parent $|\Delta E / E| \lesssim 10^{-14}$8, with MST-chain coordinates defined as

$|\Delta E / E| \lesssim 10^{-14}$9

The original coordinates are recovered by path summation to a root (com, or arbitrary particle). For evaluating vector differences $N_\mathrm{part} \sim 10^4$0, if both nodes are "close" in the tree topology (within $N_\mathrm{part} \sim 10^4$1 edges via the lowest common ancestor), their MST-path is used directly; otherwise, direct coordinates are referenced.

Key advantages of the MST decomposition relative to chain coordinates include:

• Average depth $N_\mathrm{part} \sim 10^4$2, compared to chain length $N_\mathrm{part} \sim 10^4$3; this considerably reduces floating-point accumulation errors and the arithmetic required for coordinate transforms.
• In practice, the energy drift induced by coordinate roundtrips is an order of magnitude smaller ($N_\mathrm{part} \sim 10^4$4) than in the classical chain.
• The MST structure generalizes to arbitrary $N_\mathrm{part} \sim 10^4$5 without pathological dependence on initial ordering or hierarchical clustering (Rantala et al., 2020).

3. Gragg–Bulirsch–Stoer Extrapolation and Adaptive Integration

MSTAR employs the Gragg–Bulirsch–Stoer (GBS) extrapolation method atop the regularised $N_\mathrm{part} \sim 10^4$6-time leapfrog integration, achieving rapid convergence toward the zero step-size solution. For a global timestep $N_\mathrm{part} \sim 10^4$7, a sequence of substep counts $N_\mathrm{part} \sim 10^4$8 (Deuflhard sequence $N_\mathrm{part} \sim 10^4$9) is chosen, with physical step $N$0, leading to a composition: $N$1 where $N$2 and $N$3 denote drift and kick phases. Extrapolation is realized through the Neville-Aitken triangular scheme, with error estimates as

$N$4

and convergence accepted when

$N$5

Step size $N$6 is then adapted for the next interval according to

$N$7

with $N$8. As the underlying leapfrog is symmetric, the extrapolated local error acquires $N$9 scaling, promoting extremely high conservation of energy, angular momentum, and Laplace–Runge–Lenz vector, e.g., $$H(\{ \mathbf r_i,\mathbf v_i\}) = \sum_{i=1}^N\frac{1}{2}m_i\,\mathbf v_i^2 - \sum_{i<j}\frac{G\,m_i\,m_j}{|\mathbf r_j-\mathbf r_i|},$$0, $$H(\{ \mathbf r_i,\mathbf v_i\}) = \sum_{i=1}^N\frac{1}{2}m_i\,\mathbf v_i^2 - \sum_{i<j}\frac{G\,m_i\,m_j}{|\mathbf r_j-\mathbf r_i|},$$1 radians (Rantala et al., 2020).

4. Parallelisation Architecture and Scalability

MSTAR integrates force calculation and extrapolation parallelism using the Message Passing Interface (MPI):

• Force-loop parallelism: $$H(\{ \mathbf r_i,\mathbf v_i\}) = \sum_{i=1}^N\frac{1}{2}m_i\,\mathbf v_i^2 - \sum_{i<j}\frac{G\,m_i\,m_j}{|\mathbf r_j-\mathbf r_i|},$$2 pairwise force computations are divided among $$H(\{ \mathbf r_i,\mathbf v_i\}) = \sum_{i=1}^N\frac{1}{2}m_i\,\mathbf v_i^2 - \sum_{i<j}\frac{G\,m_i\,m_j}{|\mathbf r_j-\mathbf r_i|},$$3 MPI processes. Linear strong scaling is observed up to $$H(\{ \mathbf r_i,\mathbf v_i\}) = \sum_{i=1}^N\frac{1}{2}m_i\,\mathbf v_i^2 - \sum_{i<j}\frac{G\,m_i\,m_j}{|\mathbf r_j-\mathbf r_i|},$$4, above which communication overheads dominate.
• GBS subdivision parallelism: The $$H(\{ \mathbf r_i,\mathbf v_i\}) = \sum_{i=1}^N\frac{1}{2}m_i\,\mathbf v_i^2 - \sum_{i<j}\frac{G\,m_i\,m_j}{|\mathbf r_j-\mathbf r_i|},$$5 substeps for each GBS order are partitioned among $$H(\{ \mathbf r_i,\mathbf v_i\}) = \sum_{i=1}^N\frac{1}{2}m_i\,\mathbf v_i^2 - \sum_{i<j}\frac{G\,m_i\,m_j}{|\mathbf r_j-\mathbf r_i|},$$6 groups. Each group integrates its assigned substeps independently prior to global extrapolation.
• Combined efficiency: The optimal product $$H(\{ \mathbf r_i,\mathbf v_i\}) = \sum_{i=1}^N\frac{1}{2}m_i\,\mathbf v_i^2 - \sum_{i<j}\frac{G\,m_i\,m_j}{|\mathbf r_j-\mathbf r_i|},$$7 is,

$$H(\{ \mathbf r_i,\mathbf v_i\}) = \sum_{i=1}^N\frac{1}{2}m_i\,\mathbf v_i^2 - \sum_{i<j}\frac{G\,m_i\,m_j}{|\mathbf r_j-\mathbf r_i|},$$8

which delivers near-ideal scaling; the law arises from maximizing overall speed-up before communication and synchronisation penalties set in (Rantala et al., 2020).

5. Performance, Benchmarks, and Comparison with AR-CHAIN

Metric	MSTAR	AR-CHAIN
Energy conservation	$$H(\{ \mathbf r_i,\mathbf v_i\}) = \sum_{i=1}^N\frac{1}{2}m_i\,\mathbf v_i^2 - \sum_{i<j}\frac{G\,m_i\,m_j}{|\mathbf r_j-\mathbf r_i|},$$9	$t$0
$t$1 run	Weeks on 400 CPUs (1 Gyr)	Many months
Serial speed	$t$2–$t$3 faster	—
Parallel speedup	$t$4–$t$5 (MSTAR vs. serial); up to $t$6 (AR-CHAIN baseline)	Saturates at $t$7 CPUs, then declines
Scalability	Ideal to $t$8 CPUs, saturates by $t$9	—

Serial MSTAR outperforms AR-CHAIN by $B = -H$0–$B = -H$1 due to improved cache and simpler GBS logic. Parallel MSTAR gains $B = -H$2–$B = -H$3 over its serial case and, crucially, up to $B = -H$4 over AR-CHAIN for $B = -H$5–$B = -H$6. Scalability is robust to $B = -H$7 CPUs, compared to AR-CHAIN's sub-$B = -H$8 plateau, making MSTAR suitable for simulations with unprecedented particle counts and precision (Rantala et al., 2020).

6. API-Level Coupling with Tree and Fast Multipole Integrators

MSTAR is intended as an embedded subsystem integrator within hybrid $B = -H$9-body codes using tree or fast multipole solvers (e.g., GADGET-3). In this role:

• Small, high-density "collisional" regions (typically near SMBHs) containing $s$0–$s$1 particles are passed to MSTAR for regularised integration.
• Positions and velocities are returned to the global solver at synchronised intervals.
• MSTAR’s capacity to integrate without softening allows galactic simulations with $s$2–$s$3 stars and collisional subsystems to be computed self-consistently, a key advance for capturing phenomena such as SMBH binary evolution and LISA-band gravitational wave source populations (Rantala et al., 2020).

7. Further Context: Comparison and Evolution from AR-CHAIN

The AR-CHAIN class of integrators, originally formalized by Mikkola & Merritt and subsequently extended by toolkit implementations such as SpaceHub, regularises close encounters using a linear chain of coordinate differences and combines algorithmic regularisation with compensational floating-point arithmetic and high-order symplectic steps (Wang et al., 2021). MSTAR’s MST-based coordinate system extends AR-CHAIN by reducing coordinate construction error $s$4 over the chain, removing linear-order constraints, and supporting two-level parallelism required for modern high-performance computing environments.

A plausible implication is that, for future multi-physics simulations requiring accurate collisional dynamics in multi-million particle systems, MSTAR not only removes the AR-CHAIN bottleneck but defines a new scaling regime for regularised $s$5-body computation. Energy, angular momentum, and LRL invariants track reference orbits at or near machine precision across long timescales, making the technique foundational for next-generation galactic dynamics, nuclear star clusters, and gravitational wave source studies (Rantala et al., 2020, Wang et al., 2021).