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 $dt/ds = 1/U$, eliminating the singular behaviour of $1/|\mathbf r_i-\mathbf r_j|$ as pairs approach. The transformed Hamiltonian becomes

$K = (T + B) - \ln U,$

permitting a splitting into drift (momentum only) and kick (position only) operators, on which a second-order leapfrog in $s$ is exact for Keplerian orbits. As $r_{ij} \to 0$, $dt/ds \to 0$ 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 $N-1$ edges that connect all particles with minimal total length via (e.g.) Prim's algorithm. Each edge is directed from a child $i$ to a unique parent $p(i)$, with MST-chain coordinates defined as

$$\mathbf X_i = \mathbf r_i - \mathbf r_{p(i)}, \qquad \mathbf V_i = \mathbf v_i - \mathbf v_{p(i)}.$$

The original coordinates are recovered by path summation to a root (com, or arbitrary particle). For evaluating vector differences $\mathbf r_j - \mathbf r_i$, if both nodes are "close" in the tree topology (within $N_d$ 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 $\langle L_{\rm MST}\rangle \sim O(\sqrt{N})$, compared to chain length $\sim N/2$; 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 ($\sim10\times$) than in the classical chain.
• The MST structure generalizes to arbitrary $N$ 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 $s$-time leapfrog integration, achieving rapid convergence toward the zero step-size solution. For a global timestep $H$, a sequence of substep counts $n_k$ (Deuflhard sequence $n_k = 2k$) is chosen, with physical step $h = H/n_k$, leading to a composition: $[D(h/2)\,K(h)\,D(h/2)]^{n_k},$ where $D$ and $K$ denote drift and kick phases. Extrapolation is realized through the Neville-Aitken triangular scheme, with error estimates as

$\epsilon_k \sim |T_{k,k}-T_{k,k-1}|$

and convergence accepted when

$$\frac{\epsilon_k}{|T_{k,k}| + (h/n_k)\,|T_{k,k}^\prime|} \le \text{tol}.$$

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

$$H_\mathrm{new} = a_\mathrm{GBS} \left(\frac{\text{tol}}{\epsilon_k}\right)^{1/(2k-1)} H,$$

with $a_{\rm GBS}\lesssim 1$. As the underlying leapfrog is symmetric, the extrapolated local error acquires $O(h^{2k})$ scaling, promoting extremely high conservation of energy, angular momentum, and Laplace–Runge–Lenz vector, e.g., $|\Delta E/E|\lesssim 10^{-14}$, $\theta_\mathrm{LRL} \lesssim 10^{-12}$ 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: $O(N^2)$ pairwise force computations are divided among $N_\mathrm{force}$ MPI processes. Linear strong scaling is observed up to $N_\mathrm{force}\sim 0.1N_\mathrm{part}$, above which communication overheads dominate.
• GBS subdivision parallelism: The $n_k$ substeps for each GBS order are partitioned among $N_\mathrm{div}$ groups. Each group integrates its assigned substeps independently prior to global extrapolation.
• Combined efficiency: The optimal product $$N_\mathrm{CPU} = N_\mathrm{force} \times N_\mathrm{div}$$ is,

$N_\mathrm{CPU} \sim 0.2\, N_\mathrm{part},$

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	$|\Delta E/E| \lesssim 10^{-14}$	$|\Delta E/E| \lesssim 10^{-12}$
$N_\mathrm{part}=5,000$ run	Weeks on 400 CPUs (1 Gyr)	Many months
Serial speed	$2$–$6\times$ faster	—
Parallel speedup	$55$–$145\times$ (MSTAR vs. serial); up to $1100\times$ (AR-CHAIN baseline)	Saturates at $\sim 0.02 N_\mathrm{part}$ CPUs, then declines
Scalability	Ideal to $0.1N_\mathrm{part}$ CPUs, saturates by $0.2N_\mathrm{part}$	—

Serial MSTAR outperforms AR-CHAIN by $2$–$6\times$ due to improved cache and simpler GBS logic. Parallel MSTAR gains $55$–$145\times$ over its serial case and, crucially, up to $1100\times$ over AR-CHAIN for $N\sim10^3$–$10^4$. Scalability is robust to $0.1N_\mathrm{part}$ CPUs, compared to AR-CHAIN's sub-$0.02N_\mathrm{part}$ 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 $N$-body codes using tree or fast multipole solvers (e.g., GADGET-3). In this role:

• Small, high-density "collisional" regions (typically near SMBHs) containing $N_\mathrm{reg} \sim 10^3$–$10^4$ 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 $N_\mathrm{global} \sim 10^8$–$10^9$ 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 $10\times$ 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 $N$-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).