KETJU Module: High-Precision Few-Body Dynamics

• KETJU Module is a computational framework that couples large-scale galaxy simulations with high-precision few-body dynamics around SMBHs and massive stars.
• It employs a hybrid integration scheme that divides global Tree–PM methods and local chain regularization, enabling unsoftened gravitational interactions and post-Newtonian corrections.
• The module supports parallel MPI integration and coupling with hydrodynamics, star formation, and AGN feedback, validated by benchmark comparisons with direct N-body simulations.

KETJU is a computational module designed to resolve the dynamics of point-mass objects such as supermassive black holes (SMBHs) and massive stars at sub-parsec to Schwarzschild-radius scales, within the context of large-scale galaxy and cosmological simulations. Using algorithmic regularization and direct $N$-body integration, KETJU removes the spatial resolution limitations imposed by gravitational softening in standard Tree–PM or SPH hydrodynamics codes, allowing for high-fidelity modeling of few-body interactions, post-Newtonian effects, and feedback coupling with the host galactic environment (Rantala et al., 2016, Mannerkoski et al., 2021, Johansson et al., 2022, Mannerkoski et al., 2023, Lahén et al., 2024).

1. Hybrid Architecture and Domain Decomposition

KETJU is implemented as a modular extension to codes such as GADGET-3 and GADGET-4, introducing a two-regime integration system (Mannerkoski et al., 2021, Rantala et al., 2016, Mannerkoski et al., 2023, Lahén et al., 2024):

• Global Regime: Outside of designated regions, gravity is handled by the standard Tree–PM solver (oct-tree or FMM$+$PM), coupled to SPH hydrodynamics and subgrid recipes for star formation, gas cooling, and feedback. Dynamical evolution uses softened gravity with a second-order leapfrog integrator and adaptive individual timesteps.
• Local Regularized Regions: Around each SMBH (or, in star cluster applications, each massive star and its remnant), a spherical "chain" or "KETJU region" is defined—radius $r_{\rm region} \sim 2$–$3\,\epsilon$ where $\epsilon$ is the local softening length. All relevant particles (SMBHs, stars within the region) are removed from the tree and advanced using the MSTAR (algorithmically regularized $N$-body) direct integrator, with unsoftened Newtonian or post-Newtonian (PN) gravity (Mannerkoski et al., 2021, Lahén et al., 2024).
• Coupling: At each major timestep, chain particle positions and velocities are handed off to MSTAR for direct high-precision integration in physical space, including external potential terms from the tree. After the integration, results overwrite the global code's particle data. Overlapping regions are merged to treat all close interactions simultaneously (Mannerkoski et al., 2023, Lahén et al., 2024).
• Parallelization: Chains are integrated in parallel over MPI tasks, with dynamic load balancing based on prior region CPU cost (Mannerkoski et al., 2023).

This division enables simulations to bridge length scales from $\sim 100$ kpc down to $\sim 10$ Schwarzschild radii, fully coupling SMBH/stellar microdynamics to cosmological structure formation (Mannerkoski et al., 2021, Mannerkoski et al., 2023).

2. Algorithmic Regularization and Direct Integration

KETJU uses a chain regularization technique based on the MSTAR library (an evolution of AR-CHAIN) to solve few-body interactions (Rantala et al., 2016, Lahén et al., 2024, Mannerkoski et al., 2023):

• Time Transformation: Physical time $t$ is mapped to a fictitious time $s$ using $dt/ds = \Omega(\{\mathbf{r}_i,\mathbf{p}_i\})$, e.g., $\Omega = T+B$, where $T$ is the total kinetic energy and $B$ is a function related to the system binding energy. This slows the advance of physical time during close approaches, regularizing singularities in classical $1/r^2$ gravity (Lahén et al., 2024, Mannerkoski et al., 2023).
• Chain Coordinates: Particles in each region are reordered into a minimum-spanning-tree "chain," minimizing roundoff and coordinate singularities during close encounters (Lahén et al., 2024).
• Integrator: Equations of motion in regularization time $s$ are solved using symmetric leapfrog and extrapolated to high precision via the Gragg–Bulirsch–Stoer (GBS) scheme, with per-substep tolerance $\epsilon_{\rm GBS} \sim 10^{-8}$–$10^{-7}$ and synchronization tolerance $\epsilon_t \sim 10^{-6}$–$10^{-3}$ (Lahén et al., 2024, Mannerkoski et al., 2023, Mannerkoski et al., 2021).
• Region Merging: Overlapping regularized regions are merged, ensuring that all relevant interactions (including triples and higher-order multiples) are treated in a single direct summation (Lahén et al., 2024).
• Softening Removal: All SMBH–star and SMBH–SMBH interactions are unsoftened. Optionally, a small nonzero softening can be applied for star–star encounters to stabilize the edge of the region (Mannerkoski et al., 2021, Yuan et al., 2016).
• Performance: MSTAR supports up to $\sim 10^4$ particles per region on tens of MPI ranks. For star cluster problems, integrating only the fraction of most massive stars captures most cluster physics at a fraction of the whole $N$-body computational cost (Lahén et al., 2024, Mannerkoski et al., 2021).

3. Post-Newtonian Corrections and Merger Criteria

KETJU integrates relativistic corrections up to 3.5PN order for SMBH pairs and includes PN cross-terms for $N>2$ SMBH systems (Mannerkoski et al., 2023, Mannerkoski et al., 2021):

• Pairwise PN Terms: For every SMBH pair, binary PN acceleration terms up to 3.5PN (including 1PN, 2PN, 2.5PN, 3PN, and 3.5PN) are included, covering both conservative and dissipative (gravitational-wave) effects. The 2.5PN term,

$$\mathbf{a}_{2.5\,\mathrm{PN}} = \frac{8}{5}\frac{G^2M\mu}{c^5 r^3} \dot r\, \Big[(3v^2+\tfrac{17}{3}GM/r)\hat{\mathbf{r}}-(v^2 + 3 GM/r)\tfrac{(\mathbf{r} \cdot \mathbf{v})}{r^2}\mathbf{r}\Big],$$

drives binary inspiral at sub-parsec scales (Mannerkoski et al., 2021, Mannerkoski et al., 2023).

• Three-Body PN Effects: For triplet SMBH systems, the code includes all mutual gravitational and PN-corrected forces in the direct integration. The leading-order 1PN cross-terms for $N>2$ bodies are implemented in the GADGET-4 version (Mannerkoski et al., 2023). For realistic galaxy applications, their numerical impact is validated and typically small at current resolution levels (Mannerkoski et al., 2021).
• SMBH Merger Handling: Mergers occur when the PN-predicted gravitational-wave coalescence timescale $t_c \approx -a/(4\dot a)$ falls below a threshold (e.g., $t_c < 2\,\Delta t_\mathrm{tree}$), or when the separation reaches a fixed multiple of the Schwarzschild radius (e.g., $r_\mathrm{merge} = 12R_\mathrm{S}$). The remnant’s mass, spin, and recoil kick are computed via numerical fitting formulas (Rantala et al., 2016, Mannerkoski et al., 2023).
• Lidov–Kozai Oscillations: For hierarchical triplets, KETJU resolves Kozai–Lidov cycles and their suppression by relativistic precession, with the relevant secular evolution equations included up to the necessary order (Mannerkoski et al., 2021, Johansson et al., 2022).

4. Baryonic and Stellar Physics Coupling

KETJU is embedded within the Gadget-based SPH framework (SPHGAL), allowing direct coupling to baryonic physics (Mannerkoski et al., 2021, Johansson et al., 2022, Lahén et al., 2024):

• Hydrodynamics: Standard SPH methods are used with pressure–entropy formulations and advanced kernels (e.g., Wendland C$^4$), tracking gas with metal-dependent cooling (11 elements), artificial conduction, and viscosity (Mannerkoski et al., 2021, Lahén et al., 2024).
• Star Formation and Feedback: Stochastic gas-to-star conversion applies above density thresholds (e.g., $n_\mathrm{H} = 0.1\,\mathrm{cm}^{-3}$), on local dynamical times. Feedback is modeled via energy and mass injection from supernovae, stellar winds, and AGN accretion-driven processes (Mannerkoski et al., 2021, Johansson et al., 2022).
• SMBH Seeding and Growth: Dark matter halos above a threshold (e.g., $M_\mathrm{DM} \ge 10^{10}\,h^{-1}\,M_\odot$) are seeded with SMBHs and grown via accretion models (Bondi–Hoyle–Lyttleton with Eddington cap and radiative efficiency), and merger prescriptions as detailed above (Mannerkoski et al., 2021, Johansson et al., 2022).
• Star Cluster Models: In specialized applications (e.g., SPHGAL+KETJU in dwarf galaxies), regularized integration is applied for gravitational interactions around massive stars ($m_i\gtrsim3\,M_\odot$) and their remnants, permitting accurate modeling of binary formation, core collapse, stellar ejection, and cluster dissolution without softening artifacts (Lahén et al., 2024).

5. Numerical Resolution, Performance, and Validation

• Resolution Criteria: To ensure physical reliability, especially for SMBH binary hardening, the mass ratio $m_\mathrm{BH}/m_*$ must satisfy $\gtrsim 500$–$1000$ (Mannerkoski et al., 2021, Johansson et al., 2022).
• Timestepping: The global (Tree–PM/SPH) timestep is adaptive leapfrog, while block timesteps (or GBS-controlled step selection) are used within each MSTAR region, typically $\eta_{\rm GBS} \leq 10^{-7}$ (Mannerkoski et al., 2021, Mannerkoski et al., 2023, Lahén et al., 2024).
• Scalability: For Tree–PM–dominated cases or small chains, the code scales up to $\sim 100$ cores; for chains with thousands of particles (e.g., SMBH triplets or dense clusters), nearly ideal scaling is achieved up to $\sim 1000$ cores. MSTAR's $\mathcal{O}(N^2)$ scaling in each region limits practical region size to $\sim 10^4$ particles (Mannerkoski et al., 2021, Mannerkoski et al., 2023).
• Validation: Energy conservation is maintained at $\Delta E/E \lesssim 10^{-12}$ within chain regions. Dynamical friction, binary hardening, and SMBH coalescence timescales have been benchmarked against NBODY7 and analytic models, with hardening rates $H$ and $K$ for SMBH binaries within the isolated merger-remnant ranges (e.g., $H\sim12$, $K\sim0.1$ for the inner AC binary in a cosmological simulation) (Mannerkoski et al., 2021, Rantala et al., 2016). Lidov–Kozai cycles, core–Sérsic profile core scouring, and binary ejection dynamics are robustly reproduced (Mannerkoski et al., 2023, Johansson et al., 2022).

6. Astrophysical Applications and Scientific Results

KETJU has enabled a range of high-fidelity studies in both galaxy and star cluster regimes (Mannerkoski et al., 2021, Johansson et al., 2022, Rantala et al., 2016, Lahén et al., 2024, Mannerkoski et al., 2023):

• SMBH Mergers and Triplet Dynamics: Cosmological zoom-in simulations with KETJU follow the evolution of SMBH binaries and triplets formed during multiple galaxy mergers, resolving ejections, three-body slingshot events, and hierarchical Kozai–Lidov oscillations with relativistic suppression, down to tens of Schwarzschild radii (Mannerkoski et al., 2021).
• Host Galaxy Response: SMBH binary and triplet hardening excavates stellar cores, yielding mass deficits $M_{\rm def}/M_\bullet\sim3$–$7$ and velocity anisotropy signatures consistent with core-scouring theory (Mannerkoski et al., 2023).
• Cluster Evolution: In star-by-star cluster formation models for dwarf galaxies, KETJU captures binary and higher-order collisional dynamics, producing cluster mass segregation, core expansion, and realistic disruption rates (e.g., $63\%$ of clusters disrupt in $100$ Myr) in line with direct $N$-body and observed properties; collisional corrections decrease the bound supernova fraction by $\sim1.7$ (Lahén et al., 2024).
• Coalescence Timescales: For SMBH binaries, the average post-binary-formation coalescence time is $\sim200$ Myr. Eccentricity variations within the typical range $e=0.6$–$0.95$ substantially affect merger timing due to the steep dependence in the Peters (1964) formula (Johansson et al., 2022, Rantala et al., 2016).

7. Limitations, Extensions, and Best Practices

• Current Limits: Standard KETJU modules are restricted to collisionless dynamics. Gas dynamical friction, circumbinary disk torques, and AGN feedback modules may require further code development or adoption of external recipes (Mannerkoski et al., 2023).
• Code Portability: The MSTAR library is stand-alone and can be interfaced with other grid or particle-based simulation codes (e.g., AREPO, RAMSES, ENZO) by identifying and extracting local subsystems, running regularized integration, and coupling results back (Mannerkoski et al., 2023).
• Parameter Choices: For star clusters, recommended values are $\epsilon_* \sim 0.01$–$0.1\,\mathrm{pc}$, $r_{\rm KETJU} \sim 3\,\epsilon_*$, $m_i = 3$–$8\,M_\odot$ for massive-star seeds, $\eta_{\rm GBS}=10^{-7}$, and $\eta_t =10^{-3}$ (Lahén et al., 2024).
• Energy and Error Control: Integrator settings should be chosen to ensure agreement to within $<10^{-5}$ energy error over Myr timescales, without excessive direct $N$-body step counts (Lahén et al., 2024).

In conclusion, KETJU provides a unique hybrid framework for embedding fully regularized few-body dynamics within large-scale galactic and cosmological simulations, enabling point-mass treatment and post-Newtonian accuracy in regimes previously inaccessible to grid and particle codes with gravitational softening (Mannerkoski et al., 2021, Johansson et al., 2022, Rantala et al., 2016, Lahén et al., 2024, Mannerkoski et al., 2023).