KRIOS: A new basis-expansion $N$-body code for collisional stellar dynamics

TL;DR

KRIOS presents a hybrid N-body code that blends 3D self-consistent field (SCF) relaxation with Henon-style two-body scattering to model collisional stellar dynamics without assuming spherical symmetry. The method uses Zhao1996 basis functions parameterized by $(\alpha,b)$ to flexibly fit evolving cluster potentials, and updates the 3D mean field while applying energy-conserving two-body encounters to capture diffusion and relaxation. The authors validate the approach by reproducing the radial orbit instability, core collapse in Plummer spheres (isotropic, anisotropic, and rotating), and by benchmarking timing against direct N-body codes, demonstrating effective scaling and accuracy with controlled energy and angular momentum conservation. They also demonstrate the practical advantages of 3D modeling in capturing non-spherical features, anisotropy, and rotation, with clear prospects for incorporating external tidal fields and binary dynamics in future work. The work offers a scalable, accurate framework for exploring complex globular cluster dynamics beyond spherical, orbit-averaged approximations.

Abstract

The gravitational $N$-body problem is a nearly universal problem in astrophysics which, despite its deceptive simplicity, still presents a significant computational challenge. For collisional systems such as dense star clusters, the need to resolve individual encounters between $N$ stars makes the direct summation of forces - with quadratic complexity - almost infeasible for systems with $N\gtrsim 10^6$ particles over many relaxation times. At the same time, the most common Monte Carlo $N$-body algorithm - that of Hénon - assumes the cluster to be spherically symmetric. This greatly limits the study of many important features of star clusters, including triaxiality, rotation, and the production of tidal debris. In this paper, we present a new hybrid code, KRIOS, that combines 3D collisionless relaxation using an adaptive self-consistent field method with collisional dynamics handled via Hénon's method. We demonstrate that KRIOS can accurately model the long-term evolution of clusters and provide its complete phase-space information over many relaxation times. As a test of our new code, we present detailed comparisons to well-known results from stellar dynamics: (i) the collisional evolution of a family of Plummer spheres with varying anisotropy and rotation to core collapse, and (ii) the emergence of the radial-orbit instability in radially anisotropic star clusters, including its non-spherical effects.

Figures (15)

• Figure 1: Behavior of the radial basis function, $U_n^{\ell}(r)$, of the lowest Zhao potential basis elements $n=\ell=0$, defined in Equation \ref{['eq:Unl_Zhao']}. We display four $\alpha$ values (shown in blue, red, green and purple by increasing order) and two $b$ values (0.1 in dashed lines, 1.0 in solid lines). The $\alpha$ parameter is a power-law index that sets the "cuspiness" of the potential-density profile. For example, $\alpha=1/2$ corresponds to a Plummer sphere and $\alpha=1$ to the Hernquist potential. The $b$ parameter can be qualitatively related to the scale length of the system.
• Figure 2: KRIOS uses two frames to describe the dynamics of the cluster: (i) a fixed reference frame (in black), centered at $\mathbf{0}$ and attached to a fixed set of Cartesian coordinates $(x,y,z)$; (ii) a time-dependent frame called the SCF frame (in purple), centered at the density center of the cluster $\mathbf{r}_{\mathrm{c}}(t)$, attached to a translation of the same Cartesian coordinates, and to which we add a set of spherical coordinates $(r,\vartheta,\phi)$. Frame (i) is used to integrate the dynamics in an inertial frame, while Frame (ii) is adapted to the description of the cluster and is used to expand its potential-density pair using the SCF method.
• Figure 3: Illustration of the 2-body gravitational scattering in the plane of deflection for the particles $i$ and $i+1$, which are selected by a predetermined pairing scheme (see Section \ref{['subsec:particle_pairing']}). We defined the relative position $\Delta \mathbf{r}=\mathbf{r}_i-\mathbf{r}_{i+1}$ and the relative velocity $\boldsymbol{w}=\boldsymbol{v}_i-\boldsymbol{v}_{i+1}$ of the two interacting particles. These two particles are deflected by an angle $\beta_e$ and follow hyperbolic trajectories during their interaction. This results in a modification of their velocities which exactly conserves energy, but only approximately conserves angular momentum.
• Figure 4: Schematic flowchart of the KRIOS code during one timestep. The sections describing each entries are given within brackets.
• Figure 5: Profile of the potential (top panel, $\psi$) and radial density (bottom panel, $\rho_r$) for a cluster nearing core collapse. The instantaneous potential and density are represented in green, while the SCF expansions are represented in red (using an optimized Zhao1996 family) and in purple (using the CB73 family). Although both radial expansions are limited by $n_{\max}=10$, it is clear that the CB73 expansion struggles to converge towards the instantaneous potential-density pair, whereas the optimized Zhao1996 expansion quickly converges towards it, illustrating both its necessity and effectiveness.