Kerwann Tep, Brian T. Cook, Carl L. Rodriguez, Jiya Jolly, Eddie Sawin, Michael S. Petersen, Christoph Gaffud

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\'enon - 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 a collisional dynamics part handled via H\'enon'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 two well-known results from stellar dynamics: (i) the collisional evolution of an isotropic Plummer sphere to core collapse, and (ii) the emergence of the radial-orbit instability in radially anisotropic star clusters, including its non-spherical effects.