Arepo-CD: Coupling Collisional $N$-Body Dynamics with Moving-Mesh Magnetohydrodynamics

Introduction

Star clusters form in the densest regions of giant molecular clouds (GMCs), where supersonic turbulence, self-gravity, and magnetic fields produce a highly structured interstellar medium [88, 85, 94]. Cloud collapse generates filaments, hubs, and compact overdensities in which stars form in a strongly clustered and hierarchical fashion [90, 83]. Observations identify star-forming regions at different evolutionary stages, from deeply embedded star-forming regions to exposed young massive clusters. These observations show that the transition is rapid and dynamically complex [124, 3, 123]. During this phase, stars continue to form and inject feedback as the natal cloud is dispersed [60, 13]. Cluster formation therefore depends on the complex interplay among gas dynamics, star formation, feedback, and stellar encounters over a wide range of scales. The coupled, multiscale nature of these processes makes numerical simulations essential for a self-consistent theoretical description of star cluster formation.

On the one hand, hydrodynamic simulations provide a natural framework for studying this process because the gaseous environment strongly affects the properties of forming stars [e.g. 21, 22, 79, 67, 46, 54] (H. Li et al. 2019). Current simulations can follow the formation and evolution of individual stars in molecular clouds and star-forming regions, from the initial cloud collapse to the dispersal of the natal gas (e.g. M. R. Bate 2012, 2014, 2019; C.-C. He et al. 2019; D. Guszejnov et al. [55]; M. Y. Grudić et al. [51, 52]; J. P. Farias et al. [36]). However, most such simulations use softened stellar gravity and have limited spatial resolution, which restrict their treatment of collisional stellar dynamics. This limitation artificially suppresses close encounters, binary formation and hardening, few-body scattering, dynamical ejections, and rapid mass segregation. This can alter both the internal evolution of a cluster and the locations at which massive stars inject feedback.

On the other hand, collisional $N$-body simulations provide a complementary approach to modeling stellar encounters (e.g. L. Wang et al. [141, 138]; A. Rantala et al. [119]). Here, “collisional” refers to stellar systems in which few-body gravitational encounters are dynamically important. With high force accuracy and specialized integration methods, $N$-body calculations can follow relaxation, core contraction, binary-mediated energy exchange, dynamical ejections, and long-term structural evolution (H. Baumgardt et al. [11]). They have also clarified how initially substructured stellar distributions evolve and how strong encounters reshape young clusters (S. P. Goodwin & A. P. Whitworth [49]). Their main limitation during the embedded phase is that they do not follow star formation self-consistently in an evolving gaseous environment. Conventional models begin with prescribed stellar initial conditions (ICs). They either neglect the natal gas or represent it with a static or time-dependent background potential (J. P. Farias et al. [37]; C. Schoettler et al. [125]; F. Dinnbier & P. Kroupa [31]; F. Dinnbier et al. [32]). These models therefore have to assume, rather than predict, the timing and geometry of gas dispersal.

Several hybrid “hydro + $N$-body” approaches combine hydrodynamic solvers with accurate collisional dynamics (e.g. J. E. Wall et al. [137]; B. Polak et al. [107]; M. S. Fujii et al. [43]; M. Fujii et al. 2024; Y. Jo et al. 2024; N. Lahén et al. [87]; Y. Bernard et al. [12]; E. Chung et al. [17]; Y. Jo et al. 2026). These approaches either couple specialized few-body solvers through external interfaces or embed stellar-dynamics algorithms directly within hydrodynamic codes. Interface-based methods are flexible and can combine mature solvers, but require frequent communication and synchronization between codes (e.g. J. E. Wall et al. [137]; B. Polak et al. [107]; M. S. Fujii et al. [43]; M. Fujii et al. 2024). Alternatively, embedded methods share the gravity solver, timestep hierarchy, and parallel infrastructure of the hydrodynamic code, although existing implementations differ in their treatments on collisional gravity or coupled baryonic physics (N. Lahén et al. [87]; Y. Bernard et al. [12]; Y. Jo et al. 2024; E. Chung et al. [17]; Y. Jo et al. 2026).

The main numerical challenge in such calculations is the wide range of spatial and temporal scales. Molecular-cloud evolution proceeds over several Myr, whereas tight binaries and close encounters can have characteristic timescales as short as days. Removing gravitational softening restores the force accuracy required for these systems, but produces an extremely deep timestep hierarchy in a conventional integrator. An efficient method must therefore separate compact few-body dynamics from the slower evolution of the surrounding cloud and apply high-accuracy integration only where it is needed.

Motivated by this requirement, we present Arepo-CD, a framework that couples star-by-star collisional $N$-body dynamics to the moving-mesh magnetohydrodynamics (MHD) code Arepo (V. Springel [130]; R. Pakmor et al. [102]; R. Pakmor & V. Springel [103]; R. Pakmor et al. [104]; R. Weinberger et al. [142]). The framework treats star-star gravity without softening. Interactions involving gas or collisionless components remain softened. Arepo-CD combines hierarchical time-stepping with on-the-fly identification of compact stellar groups. It evolves their internal dynamics with a fourth-order forward symplectic integrator and slow-down algorithmic regularization (S. Oshino et al. [101]; L. Wang et al. [140]; A. Rantala et al. [120, 119]). This specialized treatment is applied to dynamically demanding subsystems. The rest of the calculation remains within the native Arepo framework and is compatible with its other physics modules. We couple Arepo-CD to the Rigel radiation MHD model for star formation and feedback from individual stars (Y. Deng et al. [27, 28, 30]). This coupling allows embedded clusters to form during GMC collapse and evolve under the simultaneous influence of gas dynamics, feedback, and stellar encounters, while resolving binary formation, strong scattering, and dynamical ejections at an efficient computational cost suitable for hydrodynamic applications.

In this paper, we describe the numerical framework of Arepo-CD and present its first astrophysical applications. In Sect. 2, we describe the hybrid gravity treatment, integration scheme, timestep criteria, group-search algorithm, and coupled star formation and feedback model. In Sect. 3, we validate the method against collisional $N$-body benchmarks. In Sect. 4, we present the first GMC applications. In Sect. 5, we discuss related frameworks, numerical performance, limitations, and future directions. In Sect. 6, we summarize our conclusions.

Numerical Methods

In this section, we describe the main numerical methods implemented in Arepo-CD, which incorporates collisional $N$-body dynamics into the moving-mesh MHD code Arepo. Arepo provides an efficient tree-gravity solver, moving-mesh hydrodynamics, and a broad set of physics modules suitable for simulations of star clusters in gaseous environments (H. Li et al. 2019; Y. Deng et al. [28]; M. C. Smith [126].¹ However, the standard Arepo gravity solver is designed for collisionless $N$-body dynamics. Close binaries, hyperbolic encounters, triples, and higher-order multiples can impose extremely short dynamical timescales. Direct integration of these systems within the timestep hierarchy can therefore become inefficient or prohibitively expensive (V. Springel [129, 130]; R. Weinberger et al. [142]; V. Springel et al. [131]). Arepo-CD addresses this problem by identifying compact stellar subsystems on-the-fly and separating their internal few-body evolution from the large-scale hydrodynamics. We first outline the overall framework and hybrid gravity strategy in Sect. 2.1, after which the remainder of this section presents the main numerical ingredients of Arepo-CD.

Overview of the Framework

We begin by summarizing the modified hybrid gravity treatment, which Figure 1 illustrates. Particles are divided into collisionless and collisional components. The collisional component, including individual stars and stellar remnants, interacts internally through unsoftened Newtonian gravity, while all interactions involving the collisionless component, including gas, dark matter, and other species, use softened gravity. To implement this hybrid treatment, we modify the Barnes–Hut octree (J. Barnes & P. Hut [7]) in Arepo so that each tree node stores the collisional and collisionless mass components separately. In addition, when the pairwise dynamical timescale of two stellar particles becomes sufficiently short, they are assigned to the same collisional “group” (see Sects. 2.4.1, 2.2.4, and 2.4.2 for the timescale criterion, group definition, and group-search procedure, respectively). The internal dynamics of each group are then evolved with a specialized few-body solver (Eq. (19)). During group integration, the gravitational influence of particles outside the group, including the remaining stellar particles and collisionless components, is approximated as a constant external background. This influence is accounted for separately through the long-range kicks (Sect. 2.2.4), rather than during the group integration.

Figure 2 illustrates the Arepo-CD workflow for dynamical integration only. Below, we describe the additional steps indicated by dashed lines when the Arepo-CD module is enabled. At the end of each timestep, the force tree walk distinguishes unsoftened stellar interactions from softened interactions and also records stellar-neighbor information for subsequent group identification (Sect. 2.4.1). At the beginning of the next timestep, this information is used to identify compact stellar groups (Sect. 2.4.2). We then apply the revised timestep criteria for collisional stellar particles in the groups (Sect. 2.3). Particles in the identified groups are subsequently reassigned to the smallest timestep currently present in the timestep hierarchy to synchronize them with the other particles. The groups are then evolved in parallel using the hybrid integration scheme described in Sect. 2.2, which combines a fourth-order forward symplectic integrator and slow-down algorithmic regularization. This resolves internal group dynamics without imposing the shortest internal timesteps on the entire simulation. Finally, Sect. 2.5 describes the coupling between Arepo-CD and the Rigel star formation and feedback model used in the hydrodynamic applications.

Time Integration

Star-by-star cluster formation simulations span a wide range of spatial and temporal scales, from the global cluster evolution to close encounters and compact few-body systems such as binaries and triples. This range challenges the standard Arepo leapfrog integrator. For example, a stable tight binary with an orbital period below one day must be followed over several Myr. This requires a timestep hierarchy of order Myr/day $\sim10^{9}$. Resolving such systems directly with the default global integration framework would be inefficient, because even an update involving only a few stellar particles can require repeated force-tree walks and thereby slow the entire simulation. Moreover, the default scheme is only second-order accurate, and accumulated truncation errors over many orbits can artificially disrupt hard binaries or cause spurious cluster expansion [1]. Modern collisional stellar-dynamics codes therefore typically use fourth- or higher-order methods [141, 138, 119].

To address this problem, no single integration method can provide both high accuracy and high efficiency across the full dynamical range of collisional systems. We therefore combine complementary methods into a hybrid Hamiltonian-splitting scheme with three levels (Sect. 2.2.4): (1) leapfrog, the native Arepo scheme, advances the large-scale evolution and handles interactions between collisional stars and collisionless components (Sect. 2.2.1); (2) the forward symplectic integrator (FSI) accurately evolves the internal dynamics of compact stellar groups (Sect. 2.2.2); and (3) slow-down algorithmic regularization (SDAR) treats the stiffest close interactions efficiently, particularly hard binaries and triples (Sect. 2.2.3). To improve efficiency, independent compact stellar groups are distributed across parallel tasks and assigned group timesteps longer than their shortest internal dynamical timescales (Sect. 2.2.5).²

Leapfrog Integrator

We first briefly review the basic leapfrog integrator used in Arepo, written in operator-splitting form [129, 131, 92], as a foundation for the operator formulations introduced below. For a Hamiltonian system, time evolution over a timestep $\Delta t$ can be approximated using the standard Hamiltonian splitting $H = T + V$, where $H$ is the Hamiltonian, and $T$ and $V$ are the kinetic and potential energies, respectively [145]. The corresponding canonical maps are the drift ($D$) and kick ($K$) operators, respectively. The second-order, time-symmetric Strang splitting gives the standard kick-drift-kick (KDK) leapfrog update,

where $K(\Delta t/2)$ is a half-step “kick” that updates momenta, or equivalently velocities, under the potential force. The full-step “drift” $D(\Delta t)$ updates positions under the kinetic term. This KDK leapfrog is time-reversible and symplectic, making it efficient and robust for long-term gravitational integrations, with bounded energy oscillations rather than the secular drift often seen in non-symplectic schemes (e.g., Runge-Kutta; [129]). We note that a different Hamiltonian splitting leads to different operators, as discussed in Sect. 2.2.4, where we introduce an alternative splitting scheme.

Fourth-order Forward Symplectic Integrator

We adopt a fourth-order symplectic integrator with strictly positive substeps to improve the integration of close encounters and compact few-body systems [145, 15, 16, 143, 62, 24, 63, 120, 64, 119]. Following A. Rantala et al. [120, 119] and Y. Bernard et al. [12], we write the forward symplectic integrator (FSI) as a five-stage map over one timestep $\Delta t$, consisting of a kick, drift, gradient-corrected kick, drift, and kick. Using the drift operator $D(\Delta t)$, the kick operator $K(\Delta t)$, and the gradient-corrected kick operator $\widetilde{K}(\Delta t)$, this cycle can be written as

The kick, drift, and gradient-corrected kick update the particles as follows:

Here $a_i$ is the Newtonian gravitational acceleration of particle $i$ and $G$ is the gravitational constant. We define $r_{ji}\equiv r_j-r_i$ with $r_{ji}\equiv\lvert r_{ji}\rvert$, and $a_{ji}\equiv a_j-a_i$.

The gradient-corrected kick in Eq. (5) makes the FSI composition fourth-order accurate. Its $O(\Delta t^3)$ correction cancels the leading error terms of the second-order scheme. Compared with high-order symplectic integrators in H. Yoshida [145], the fourth-order FSI uses strictly positive sub-step lengths, which is advantageous for irreversible dynamics (e.g., with gravitational-wave dissipation) and can improve the efficiency of hierarchical recursive integration. The method also performs better for the Kepler problem than schemes with negative substeps (S. A. Chin [15]).

Slow-Down Algorithmic Regularization

Increasing the integration order does not remove all difficulties associated with close encounters and compact few-body systems, since highly eccentric binary stars and close hyperbolic encounters can still impose prohibitively small timesteps. We therefore include slow-down algorithmic regularization³ (SDAR; L. Wang et al. [140]), a few-body integration method that combines two complementary methods. The slow-down (SD; S. Mikkola & S. J. Aarseth [96]) technique reduces the computational cost for weakly perturbed hard binaries by artificially slowing their internal orbital phase while preserving the secular effects of external perturbations. Algorithmic regularization (AR; S. Mikkola & K. Tanikawa [97]; M. Preto & S. Tremaine [115]) is a time-transformed explicit symplectic scheme that improves robustness and conservation during close encounters.

For a perturbed binary, the SD technique rescales only the internal binary motion while leaving the center-of-mass (CoM) motion unchanged. In the laboratory frame, the equations of motion for component $i=1,2$ with companion $j\ne i$ are

where $v_c$ is the CoM velocity, $V_{\mathrm{ext}}(r)$ is the external potential energy, and $\kappa$ is the SD factor. This factor $\kappa$ slows only the internal motion in the CoM frame, effectively averaging external perturbations over multiple physical orbits and allowing larger timesteps. It is set from the relative strengths of the internal binary interaction and external perturbations (L. Wang et al. [140]):

with $C$ a scaling constant, $a$ and $e$ the binary semi-major axis and eccentricity, $N_p$ the number of perturbers, and $r_c$ the CoM position. We also require the slow-down period $T_{\mathrm{SD}}=\kappa T$ to satisfy $T_{\mathrm{SD}}\le\Delta t$.

Algorithmic regularization is formulated by extending phase space and applying a time transformation. We treat physical time $t$ as an additional coordinate, with conjugate momentum $p_t=-H$. The extended coordinate, momentum, and phase-space vectors are $R\equiv(r,t)$, $P\equiv(p,p_t)$, and $W\equiv(R,P)$, respectively. We also define the extended kinetic and potential energies as $T(P)\equiv T(p)+p_t$ and $V(R)\equiv V(r,t)$. AR introduces a fictitious time $s$ via the time-transformation function $g$:

This choice gives an extended Hamiltonian written in separable logarithmic form,

which permits explicit symplectic integration, for example with leapfrog (Eq. (1)). The corresponding equations of motion are

Together, SD and AR efficiently evolve embedded few-body subsystems within larger $N$-body simulations while maintaining accuracy over long times.

Hamiltonian Splitting through the Particle-Particle Particle-Tree (P$^{3}$T) Method

We combine the above integrators into a hybrid Hamiltonian-splitting scheme to achieve both high accuracy and efficiency across the full dynamical range of collisional systems. Unlike the standard splitting $H = T + V$, this approach uses the separation of dynamical timescales in collisional $N$-body systems. Most particles evolve smoothly on longer timescales and can be advanced efficiently with lower-order integrators, while a small subset forms compact subsystems with rapid motion requiring higher-order or regularized integration.

More explicitly, we consider an $N$-body system containing $N_{\star}$ collisional stellar particles and $N-N_{\star}$ collisionless particles. At each synchronization time, we partition the particle index set $\{1,2,\ldots,N\}$ into disjoint subsets. The subsets $\mathcal{G}_k$ ($k = 1,2,\ldots,N_G$) contain stellar particles assigned to compact subsystems, which we call “groups”, and the subset $S$ contains the remaining stellar particles and all collisionless particles, which we call “singles”. This partition gives

The group-identification procedure is described in Sect. 2.4.2.

Following the P$^{3}$T-type hybrid method ([101, 138]), Arepo-CD uses this partition to split the Hamiltonian into short- and long-range components. At fixed group membership, this decomposition is algebraically exact:

with

Here, $\phi(r_{ij},\epsilon_{ij})$ denotes a general softened inverse-distance function with an effective pairwise softening length $\epsilon_{ij}$, recovering $1/r_{ij}$ as $\epsilon_{ij} \to0$. The pairwise changeover function $W_{ij}$ assigns each interaction to the short- or long-range Hamiltonian. Arepo-CD adopts a sharp changeover based on group membership,⁴ defined as

The short-range Hamiltonian $H_S$ contains the kinetic energy of all particles and the stellar gravitational interactions within each group, while the long-range Hamiltonian $H_L$ contains all remaining stellar interactions and all gravitational interactions involving collisionless particles. Following the operator notation of Eqs. (1) and (2), we denote their corresponding canonical maps over a timestep $\Delta t$ by $E_S(\Delta t)$ and $E_L(\Delta t)$, respectively. For the disjoint groups $\{\mathcal{G}_k\}_{k=1}^{N_G}$, the short-range map factorizes into independent group integrations and a drift of the single particles:⁵

where $D_S$ drifts the singles and $E^{(\mathcal{G}_k)}$ advances the internal dynamics of group $\mathcal{G}_k$. The long-range Hamiltonian $H_L$ accounts separately for the external gravitational field. Since $H_L$ contains only potential terms, its map is a pure kick, $E_L = K_L$.

We hierarchically couple the three integration methods introduced above through two levels of the short- and long-range Hamiltonian splitting. (1) At the outer level, second-order KDK leapfrog integrates collisionless particles and slowly evolving collisional stellar particles (Sect. 2.2.1), and fourth-order FSI evolves the identified stellar groups (Sect. 2.2.2). (2) At the inner level, fourth-order FSI evolves most particles within each group, and SDAR integrates the identified tight subgroups (Sect. 2.2.3). Further details are given below.

Using a symmetric Strang composition, analogous to Eq. (1), the outer coupling gives the full evolution over one timestep $\Delta t$ as

Interactions between stellar groups and single particles outside them are handled by the long-range kicks $K_{\mathrm{L}}$, evaluated using the gravitational force octree (Sect. 2.1; V. Springel [129, 130]; R. Weinberger et al. [142]; V. Springel et al. [131]). Single particles outside the groups retain the standard KDK evolution without group integration.

To balance accuracy and computational efficiency within each stellar group, we identify tight subgroups, such as close binaries and hierarchical triples, using a partition similar to Eq. (11). Fourth-order FSI integrates the general few-body dynamics, while SDAR handles the stiff subgroups more efficiently. For compactness, we define the combined SDAR evolution of all tight subgroups within group $\mathcal{G}_k$ as

where $\{\mathcal{F}_l\}_{l=1}^{N_{\mathrm{F}}^{(\mathcal{G}_k)}}$ denotes the set of tight subgroups identified within $\mathcal{G}_k$, and the map $E_{\mathrm{SDAR}}^{(\mathcal{F}_l)}(\Delta t)$ advances subgroup $\mathcal{F}_l$ with SDAR over one timestep $\Delta t$. The full evolution map of group $\mathcal{G}_k$ over $\Delta t$ is then

where $K_{\mathrm{L}}^{(\mathcal{G}_k)}$ and $D_{\mathrm{S}}^{(\mathcal{G}_k)}$ denote the kick and drift operators, respectively, acting on the remaining single stellar particles within group $\mathcal{G}_k$. The operator $\widetilde{K}_{\mathrm{L}}^{(\mathcal{G}_k)}$ is the corresponding gradient-corrected kick. For sufficiently small groups, evolving the entire group with SDAR is more efficient. We therefore use SDAR alone when $|\mathcal{G}_k| \leq N_{\mathrm{SDAR,max}}$. Following L. Wang et al. [140], we typically set $N_{\mathrm{SDAR,max}} \leq4$.

Substituting the compact stellar group integration of Eq. (19) into the global time evolution of Eq. (17) yields the three-level evolution operator that forms the core integration algorithm of Arepo-CD. The combined scheme is formally second-order accurate because the outermost coupling uses a second-order KDK leapfrog integrator. Its numerical accuracy is also affected by the approximate gravitational forces from the octree. Nevertheless, the few-body dynamics are treated with higher-order integrators, while the three-level structure covers a wide dynamical range. The benchmark tests in Sect. 3 show that Arepo-CD agrees well with $N$-body calculations and reproduces the key dynamical properties of star clusters.

Parallel Group Integration

Next, we provide further details on the integration of the compact stellar groups $\{\mathcal{G}_k\}_{k=1}^{N_{\mathrm{G}}}$ defined in Sect. 2.2.4. For each group, we compute its CoM properties and transform its members’ positions and velocities into the CoM frame. This separates the CoM motion from the internal motion. The CoM motion, together with that of single stellar and collisionless particles, is advanced with the standard KDK leapfrog scheme (Eq. (1)), while the internal relative dynamics are evolved with the fourth-order FSI+SDAR solver (Eq. (19)). The global timestep is determined from active particles outside the groups together with the group CoM contributions. All members of a group are assigned the same global timestep as their CoM to maintain synchronization with the global KDK integration. The above approach is also applied to subgroups integrated with SDAR within each compact stellar group.

The decomposition into disjoint group integrations in Eq. (17) also allows independent group integration in parallel across MPI ranks, as external gravitational interactions are handled separately through the long-range kick $K_{\mathrm{L}}$ rather than during internal group integration. Arepo-CD distributes the groups with a greedy load-balancing scheme based on an approximate integration-cost estimate following L. Wang et al. [140]. Because only group-member data are exchanged during this stage, the communication overhead between MPI ranks remains minimal. After integration, the updated particle properties are returned to the MPI ranks on which the corresponding particles reside.

As an accuracy check, the total internal energy $E$ of each group is recomputed at the end of the integration. If the relative energy error exceeds a predefined threshold, the ICs of the group are written to a debug file for further inspection. An optional configuration also records the cumulative integration error over all groups for performance diagnostics and parameter tuning.

Timestep Criteria

We introduce timestep criteria for collisional stellar particles to resolve short dynamical timescales accurately (W. Dehnen & J. I. Read [25]) and to support the group identification in Sect. 2.4. The standard Arepo gas timestep criteria already account for the relevant hydrodynamic and MHD constraints and are therefore left unchanged (V. Springel [130]; R. Weinberger et al. [142]).

For stellar particles, however, the default gravitational timestep is only controlled by the acceleration criterion commonly used in cosmological simulations (C. Power et al. [114]),

Here $\epsilon$ is a characteristic length scale, which we set to the gravitational softening length. The particle acceleration is $\mathbf{a}$ and $\eta_{\mathrm{accel}}$ is a dimensionless safety parameter. Although this criterion is effective in many collisionless applications, it is not sufficiently robust for the collisional stellar dynamics followed in Arepo-CD. Close binaries, higher-order multiples, and strong encounters can change stellar velocities on very short timescales [132], which are not always captured by an acceleration criterion alone. Moreover, geometric cancellation can also reduce the magnitude of the instantaneous acceleration $a$, which permits excessively large timesteps in dynamically important regions (see Figure 1 of M. Y. Grudić & P. F. Hopkins [53]). Robust and efficient timestep control in collisional star-gas simulations therefore requires additional conservative criteria for stellar particles. We introduce these criteria below.

Tidal Criterion

We first add the tidal timestep criterion of M. Y. Grudić & P. F. Hopkins [53] to avoid failures caused by geometric cancellation of the acceleration. This criterion estimates the local dynamical timescale from the tidal tensor $T = \nabla a$, i.e., the spatial gradient of the gravitational acceleration. The invariant Frobenius norm of $T$ gives

which recovers the usual Keplerian dynamical timescale in the appropriate limit. Here $\eta_{\mathrm{tidal}}$ is a dimensionless safety factor. When a softened force law is used, the tidal tensor must be computed consistently from the softened gravitational potential (P. F. Hopkins [66]; M. Y. Grudić et al. [51]).

Free-fall and Fly-by Criteria

To capture the short timescales associated with sudden velocity changes during close few-body interactions, we supplement the timestep control with two pairwise criteria based on physically motivated two-body timescales (A. Rantala et al. [120, 119]). For each stellar particle $i$, we consider neighboring stellar particles $j$ within the search radius $r_{\mathrm{search}}$, an adjustable input parameter.⁶ The two criteria use the pairwise free-fall timescale $\Delta t_{\mathrm{ff},ij}$ and fly-by timescale $\Delta t_{\mathrm{fb},ij}$, respectively. We take the minimum over neighboring stars within $r_{\mathrm{search}}$ and impose

and

Here $r_{ij}=|\mathbf{r}_{ij}|=|\mathbf{r}_i-\mathbf{r}_j|$ is the pairwise separation, $v_{ij}=|\mathbf{v}_{ij}|=|\mathbf{v}_i-\mathbf{v}_j|$ is the relative speed, and $m_{ij}=m_i+m_j$ is the total mass. The dimensionless safety factors $\eta_{\mathrm{ff}}$ and $\eta_{\mathrm{fb}}$ are specified as input parameters. In practice, we recommend $\eta_{\mathrm{ff}}=\eta_{\mathrm{fb}}\simeq2\eta_{\mathrm{tidal}}\simeq0.02$, so that the two-body criteria then primarily control the timestep, and that the tidal criterion provides an additional safeguard.

We note that, to resolve close stellar neighbors, all tree nodes that overlap the search region of radius $r_{\mathrm{search}}$ are opened. We also reuse this neighbor information for timestep calculations to identify groups during the gravity octree walk (Sect. 2.4.1).

Group and Subgroup Identification

We now describe how AREPO-CD identifies compact stellar “groups” (defined by Eq. (11)) for integration with the few-body solvers in Sect. 2.2.4. At the beginning of each timestep, the algorithm (1) identifies close stellar “pairs” (defined below) during the gravity calculation using the Barnes–Hut octree and then (2) merges these pairs hierarchically and adaptively into disjoint groups. The same procedure is used to identify subgroups within each group, except that the Barnes–Hut octree is replaced by direct force calculations. The following subsections therefore describe both stages using group identification as an example, with subgroup-specific details given where necessary.

Pair Identification

During the gravity octree walk, the pairwise free-fall and fly-by timescales, $\Delta t_{\mathrm{ff},ij}$ and $\Delta t_{\mathrm{fb},ij}$ (defined in Eqs. (22) and (23)), are reused for group identification. Two stellar particles ($i,j$) are identified as a collisional “pair” if either characteristic interaction timescale is sufficiently short. Specifically, we require

or

Here $r_{ij}$, $v_{ij}$, and $m_i+m_j$ have the same definitions as in Eqs. (22) and (23) and we require $j\ne i$ and $r_{ij}<r_{\mathrm{search}}$. The input parameter $\Delta t_{\mathrm{group}}$ identifies stellar particles with short dynamical timescales that require specialized few-body integration, with a typical value of $\sim10\ \mathrm{yr}$.

The pair identification is performed during the timestep calculation in the gravity force-tree walk, avoiding an additional stellar neighbor search. Moreover, as discussed in Sect. 2.2.5, stellar particles assigned to groups no longer directly impose their short internal timescales on the global integration. The few-body solvers evolve these particles on their short internal timescales, allowing a larger global timestep and improving the overall integration efficiency.

Adaptive Group Construction

As Figure 3 illustrates, the identified collisional pairs are used to construct groups in a friends-of-friends (FOF) way [65, 119, 12]. The two stars in each pair are assigned to the same collisional group. If either star forms another pair with a third star, the third star is added to the same group. This procedure is implemented using a disjoint-set data structure and the union–find algorithm, and is repeated until no additional stars can be added. The group search first identifies groups for the fourth-order FSI, and then identifies smaller subgroups within each FSI group that require SDAR integration. The first stage is parallelized across MPI tasks, whereas the second is performed locally because all members of an FSI group are already stored on the same CPU. For the parallel search, we distinguish between “local pairs”, whose two members reside on the same MPI task, and “remote pairs”, whose members reside on different tasks. Each task first constructs groups from its local pairs. The tasks then exchange remote pairs and use all pairs to merge groups across task boundaries, yielding the final disjoint groups and ungrouped single stars.

This approach provides an alternative to the default FOF algorithm in Arepo with lower overhead [130, 142]. The standard Arepo FOF method requires an additional tree containing only stars and a separate tree traversal, introducing extra memory and computational costs. By contrast, our approach identifies stellar pairs directly during the gravity calculation. The pairs must be stored explicitly, but the additional memory cost is modest because the collisional groups considered here are much smaller than the halo-scale FOF groups commonly used in cosmological simulations (e.g., [23]). For example, even if all stellar particles belong to isolated binary stars, the memory cost remains $O(N)$. In practice, the group-search cost is negligible compared with the subsequent group integration. For an extremely compact cluster, the memory required to store pairs can approach $O(N^2)$, although this can be alleviated by the adaptive group-search strategy described below.

More importantly, such rare but prohibitively expensive groups also affect the group integration, since the largest group usually dominates the integration runtime [141, 138]. The strategy is to reduce $\Delta t_{\mathrm{group}}$ to form smaller groups, together with an explicit upper limit on the size of each collisional group. Exceptionally large groups can occasionally exceed $\sim10^3$ members for the fourth-order FSI or $\sim10$ members for SDAR, making the integrators inefficient. During the union–find construction, we therefore track the size of each group and process pairs in increasing order of their minimum dynamical timescale, $\min\{\Delta t_{\mathrm{ff},ij}, \Delta t_{\mathrm{fb},ij}\}$, with $j \ne i$ and $r_{ij} < r_{\mathrm{search}}$. A merger is allowed only if the resulting group or subgroup remains below the user-defined limits, $N_{\mathrm{FSI,max}}$ (typically $\sim100$) for FSI groups or $N_{\mathrm{SDAR,max}}$ for the SDAR subgroups (introduced in Sect. 2.2.4). If a merger would exceed the group-size limit, it is skipped and the remaining candidate pairs are processed. Interactions excluded from group integration remain subject to the timestep criteria in Sect. 2.3, which can then reduce the global timestep accordingly. This trades smaller timesteps for a predictable computational cost of the embedded (sub)group integration using the (FSI+)SDAR solver.

Coupling to Rigel

To model individual star formation and feedback during star cluster formation and early evolution in a realistic gaseous environment, we couple Arepo-CD to Rigel [27, 30], a star-by-star galaxy formation framework built on Arepo. The gas is evolved with the moving-mesh hydrodynamics and ideal-MHD solvers of Arepo, while radiation is coupled through the moment-based Arepo-RT solver [146] (R. Kannam et al. 2019). The model includes detailed interstellar medium (ISM) cooling, Hii-region feedback, and supernova (SN) feedback ([78, 26, 77]). Stellar lifetimes, photon production rates, mass-loss rates, and wind velocities are determined from the initial mass and metallicity of each star using the stellar-model library described in Y. Deng et al. [27]. An updated individual-star formation model (Y. Deng et al. [30]) samples stars over the full stellar initial mass function (IMF) from 0.08 to 150 $M_\odot$ and reproduces the observed maximum stellar mass–cluster mass relation. Further details of the Rigel framework are given in Y. Deng et al. [27, 28, 30].

Arepo-CD and Rigel are coupled flexibly, allowing both the Rigel model and collisional dynamics in Arepo to be switched on or off independently for the tests presented in Sects. 3 and 4. Without Arepo-CD, stellar gravity is computed with the Barnes–Hut octree ([7]) and gravitational softening ([130, 142]). Enabling collisional dynamics also involves a modification to the initial conditions of newly formed stars. The original Rigel model randomly offsets the initial position of each newly formed star from that of its parent reservoir particle, which represents the protostellar phase through a subgrid model (see Y. Deng et al. [30] for details). We additionally introduce a random velocity perturbation to reduce the formation of unphysically compact primordial binaries. For each Cartesian coordinate, the velocity offset is drawn independently from a zero-mean Gaussian distribution with the Box–Muller method. The newly introduced velocity dispersion $\sigma_v$ is set to $\sim0.1$ km/s in our fiducial run to minimize its effect on the final stellar velocities.

Although the Rigel framework models individual stars, resolving full collisional dynamics for every star can be prohibitively expensive in galaxy-scale simulations ([87 a]). We therefore provide an option that enables collisional dynamics only in regions around massive stars (e.g. $> 8\,M_\odot$, with the threshold treated as a free parameter). This option reduces the computational cost and focuses the collisional calculation on massive stars. We disable this option in the production runs presented here because these runs cover only GMC and star cluster scales.

Benchmark Tests

We perform a suite of controlled $N$-body benchmark tests to assess the numerical methods described in Sect. 2. We compare Arepo-CD with the collisional $N$-body code PeTar⁷, a code developed for collisional stellar systems ([70, 100, 138]). We use two complementary tests. The first compares the collisional evolution with PeTar and theoretical expectations, while the second isolates the effects of collisional dynamics through comparison with simulations that use softened gravity. Throughout this work, simulations performed with PeTar are referred to as “PeTar”, those with Arepo including collisional dynamics as “Arepo-CD”, and those with softened gravity as “Arepo”.

The initial conditions (ICs) are isolated, spherically symmetric star cluster models generated with a new version of the McLuster code⁸ ([84, 139]). We adopt the H. C. Plummer [106] sphere for the initial density profile and initialize the systems in virial equilibrium. For simplicity, we neglect external tidal fields, stellar evolution, and primordial binaries. Sects. 3.1 and 3.2 give further details of the two simulation setups.

Core and Lagrangian Radii

We first compare Arepo-CD with PeTar using an idealized equal-mass star cluster whose theoretical evolution is well established. We use a cluster with $N = 10^3$ stars, each with a mass of $m_\star= 8\ M_\odot$.⁹ The initial half-mass radius is $r_{\mathrm{h},0} = 1.0$ pc. To quantify stochastic variation, we generate eight random realizations of the ICs for each code.

For this model, the initial half-mass relaxation time is estimated as

where we adopt the conventional Coulomb logarithm $\ln\Lambda= \ln(0.4N)$ (L. Spitzer [128]). For an isolated equal-mass Plummer model, core collapse occurs after $\sim14$–$17t_{\mathrm{hr}}$ (H. Cohn [18]; M. Giersz & D. C. Heggie [47]; J. Makino [91]), corresponding to $t_{\mathrm{cc}} \sim53$–$65$ Myr here. Hard-binary formation then halts core collapse and drives post-collapse expansion.

The half-mass radius approaches the scaling $r_{\mathrm{h}} \propto t^{2/3}$ (M. Giersz & D. C. Heggie [47]; M. Hénon [61]).

We now compare these expectations with the numerical evolution. Figure 4 compares the core radius and the Lagrangian radii enclosing 25%, 50%, and 75% of the bound stellar mass. The core radius and 25% Lagrangian radius probe the central collisional dynamics, while the 50% and 75% radii trace the global mass distribution and secular expansion. Before core collapse, the core contracts while the half-mass radius remains nearly constant as expected from two-body relaxation (M. Giersz & D. C. Heggie [47]). Arepo-CD is statistically consistent with PeTar throughout the evolution. Both show the expected core collapse at $t_{\mathrm{cc}} \sim53$ Myr and subsequent self-similar expansion, $r_{\mathrm{h}} \propto t^{2/3}$. The agreement across all radii shows that Arepo-CD reproduces both the inner collisional evolution and the global relaxation behavior, with residual differences smaller than the variation between realizations.

cussed in Sect. 4.4. The softened Arepo runs produce fewer high-velocity escapers, particularly in the more massive cluster. In contrast, Arepo-CD and PeTar show nearly identical speed distributions, including the high-velocity tail generated by strong collisional interactions.

Radial Density Profile and Velocity Distribution

For the second test, we consider isolated young star clusters with a realistic stellar IMF to isolate the effects of collisional dynamics. Individual stellar masses are drawn from the P. Kroupa [82] IMF, with the maximum stellar mass set from the cluster mass following Z. Yan et al. [144]. See also N. Lahén et al. [87] and Y. Deng et al. [30]. The initial three-dimensional half-mass radius is assigned from the observed cluster mass-radius relation (M. Marks & P. Kroupa [93]):

Here we adopt $\beta= 0.18 \pm0.028$ and $R_4 = 2.365 \pm0.106$ from G. Brown & O. Y. Gnedin [14], together with $f_{\mathrm{h}} \sim1/6$ to represent the compact birth radii inferred for embedded clusters observationally (e.g., A. Adamo et al. [2]). Recent dynamical models also adopt such compact radii (e.g., A. Rantala et al. [117, 116, 118]). We study two models containing $\sim10^3$ and $\sim10^4$ stars, with total stellar masses of $\sim600$ and $\sim6000\,M_\odot$, half-mass radii of 0.18 and 0.28 pc, and maximum stellar masses of $\sim15$ and $\sim55\,M_\odot$, respectively. These compact clusters have short crossing and relaxation times and are therefore expected to undergo substantial structural evolution within $\sim10$ Myr. For each model, we generate 10 random realizations.

Figure 5 compares the stellar radial density profiles and velocity distributions at $t = 10$ Myr. The top two panels show the radial density profiles of the two cluster models. The Arepo runs with softened stellar gravity remain close to the ICs. They retain a centrally concentrated stellar distribution and develop only a weak extended tail. In contrast, both Arepo-CD and PeTar deviate substantially from the ICs and the Arepo results with softened gravity, showing reduced central densities and extended outer density tails. Arepo-CD agrees well with PeTar over most of the resolved radial range. Larger stochastic differences occur mainly in the innermost bins, where only a small number of stars are present and strong few-body encounters dominate. This evolution is consistent with collisional dynamics. Two-body relaxation and mass segregation drive energy exchange in the core, while hard-binary and few-body interactions heat the cluster and move stars to larger radii (D. Heggie & P. Hut [59]).

Next, we compare the velocity distributions in the bottom two panels of Figure 5. We consider only stars with no neighbor within 0.01 pc. This selection excludes large instantaneous velocities caused by internal binary motion and focuses on escapers. We classify stars with 10 km/s $< v \leq$ 30 km/s as walkaways and those with $v > 30$ km/s as runaways, with these ranges marked by gray shaded regions ([35, 125]). Their origin is discussed in Sect. 4.4.

Self-consistent simulations of star cluster formation and dynamical evolution

Next, we turn to the full hydrodynamical GMC simulations that follow the formation and early evolution of star clusters. We analyze a fiducial model in which an isolated, isothermal, massive, turbulent GMC collapses to form stars, allowing us to characterize both the new $N$-body scheme and the resulting stellar system. Throughout this section, “RIGEL” denotes simulations with the RIGEL model described in Sect. 2.5, while “Arepo-CD” denotes simulations with RIGEL and collisional dynamics.

Initial Conditions and Simulation Parameters

The simulations start from an isolated, isothermal, turbulent GMC generated with MakeCloud¹⁰ (M. Y. Grudić & D. Guszejnov [50]), following a setup similar to Y. Bernard et al. [12]. The cloud has mass $M_{\mathrm{GMC}} = 10^5\,M_\odot$ and radius $R_{\mathrm{GMC}} = 8.8$ pc, embedded in an ambient medium with density $10^{-3}$ times that of the cloud. We place the cloud at the center of a periodic box with side length $\sim200$ pc. The gas is initialized at $T = 20$ K with mean molecular weight $\mu= 2.35$. We add a Gaussian random velocity field with power spectrum $P_k \propto k^{-2}$, scaled to a turbulent virial ratio $\alpha_{\mathrm{turb}} = 2$. We use a gas mass resolution of $m_{\mathrm{gas}} \simeq1\,M_\odot$, and follow the cloud through the embedded phase toward gas removal up to $t_{\max} \simeq2.51\,t_{\mathrm{ff}} \simeq3.42$ Myr. We define the characteristic free-fall time as

which differs slightly from that of M. Y. Grudić et al. [54]. We omit external driving and initialize a uniform magnetic field whose strength is set following T. C. Mouschovias & L. Spitzer [99]:

where $E_G$ and $E_B$ are the gravitational and magnetic energies of the gas within the cloud, respectively. We adopt $\mu_0 = 4.2$, which corresponds to $\sim21\,\mu\mathrm{G}$. Moreover, we initialize $r_{\mathrm{search}}$ to 0.01 pc $\sim2 \times10^3$ AU to better capture wide binaries (G. Duchêne & A. Kraus [34]) and close encounters.

Early Spatial Evolution

We first examine how collisional dynamics affects the early spatial evolution of the cloud and the star cluster that forms within it. Figure 6 compares the gas surface density, $\Sigma_{\mathrm{gas}}$, and mass-weighted gas temperature, $T_{\mathrm{gas}}$, in the Arepo-CD and Rigel runs. The maps show $t = 1.37$ Myr $\simeq1.0\,t_{\mathrm{ff}}$ and $t = 2.06$ Myr $\simeq1.5\,t_{\mathrm{ff}}$, together with gray stellar projections overlaid. At both epochs, the simulations show broadly similar global evolution. Gas filaments feed star formation in dense regions and hot, low-density cavities emerge as stellar feedback becomes important. At $1.0\,t_{\mathrm{ff}}$, the gas and stellar distributions are nearly indistinguishable between the two runs. By $1.5\,t_{\mathrm{ff}}$, modest local differences appear in the stellar distribution, while the overall cluster morphology remains similar. These differences are more apparent in the temperature maps. Stellar dynamics redistributes individual feedback sources, and the locations and shapes of heated regions begin to differ. Thus, collisional dynamics does not strongly alter the large-scale cloud evolution at these early times, but introduces local differences in the stellar distribution and the resulting spatial pattern of feedback (see also [12, 87]).

Star Formation Efficiency and Star Formation Rate

We next quantify the global star formation histories (SFHs) of the two runs through the star formation rate (SFR) and star formation efficiency (SFE). Here we compute the instantaneous SFR from the discrete simulation outputs and divide the stellar mass formed between two successive snapshots by

the time interval. The integrated SFE $\epsilon_\star$, which tracks the evolutionary stage of the cloud, is defined as

where $M_{\mathrm{gas}}$ is the total gas mass of the GMC. Figure 7 shows the SFHs of the two runs, following a presentation similar to that of H. Li et al. (2019). The top panel shows the same qualitative evolution in both runs. After an initial delay, the SFR rises rapidly, reaches a broad maximum of order $10^{-2}$–$10^{-1}\ M_\odot/\mathrm{yr}$ around $t \sim1$–$1.5\ t_{\mathrm{ff}}$. It then declines as stellar feedback disperses the dense gas. After $t \sim1.5\ t_{\mathrm{ff}}$, the SFR decreases steadily to $\lesssim10^{-3}\ M_\odot/\mathrm{yr}$, with increasing short-timescale variability. The integrated SFE $\epsilon_\star$ (bottom panel) rises rapidly during the main star-forming phase. It reaches $\sim0.2$ near the onset of the SFR decline, and then gradually approaches an asymptotic value of $\sim0.26$. Overall, Arepo-CD and Rigel show very similar SFR and SFE evolution, with only modest differences around the peak SFR and during the late-time decline. This agreement indicates that large-scale cloud collapse and stellar feedback primarily regulate the global SFH, while the different treatments of stellar dynamics have only a secondary effect on the integrated star formation properties.

Velocity Distribution and Runaway Stars

We next compare the stellar kinematics in the collisional Arepo-CD and collisionless Rigel runs at $t = 2.61\ \mathrm{Myr}$ ($1.91\ t_{\mathrm{ff}}$). Figure 8 shows the speed distributions for isolated stars and the massive-star subsample ($M_\star> 8\ M_\odot$). We include only isolated stars with no neighbor within $0.01\ \mathrm{pc}$ as in the bottom panels of Figure 5. The bulk distributions are similar and peak at a few $\mathrm{km/s}$, indicating comparable large-scale cluster kinematics. The main difference is the high-velocity tail. Arepo-CD produces numerous fast escapers with speeds reaching $\gtrsim100\ \mathrm{km/s}$, whereas the Rigel distribution declines sharply above $\sim10$–$20\ \mathrm{km/s}$. Adopting the conventional threshold of $v \gtrsim30\ \mathrm{km/s}$ for runaway stars ([35, 125]), the collisional treatment therefore produces a substantially larger runaway population. This is consistent with close few-body encounters that transfer binding energy into stellar kinetic energy ([128]).

Figure 9 shows the spatial distribution of these escapers at $t \simeq2.06\ \mathrm{Myr}$ ($\simeq1.5\ t_{\mathrm{ff}}$). Velocity vectors for stars with $v > 30\ \mathrm{km/s}$ are overlaid on the projected mass-weighted gas temperature, $T_{\mathrm{gas}}$. Arepo-CD produces numerous high-velocity stars that move outward from the central star-forming region, whereas Rigel has no comparable population at the same epoch. Because these systems are too young for SN disruption of binaries, these runaways arise from the dynamical ejection processes, including close few-body encounters and subcluster interactions ([113, 41, 37, 108]).

The massive-star subsamples shown in Figure 8 and Figure 9 further demonstrate that the highest velocities are not confined to low-mass stars. This is important because the ejection of massive stars can redistribute stellar feedback. Without such ejections, ionizing sources remain more deeply embedded in the molecular cloud (right panel of Figure 9), whereas runaway massive stars can carry ionizing radiation, winds, and later SN feedback into lower-density environments (left panel of Figure 9). On galactic scales, this may alter the spatial distribution and escape of ionizing radiation, while displaced massive stars can inject feedback into lower-density gas. There, feedback can couple more efficiently to the ambient ISM ([5, 6]).

Discussion

Comparison with Existing Hybrid Methods

Coupling the collisional $N$-body dynamics of individual stars to hydrodynamic star formation simulations remains technically challenging. Only a few studies address this problem. We compare the performance and results of Arepo-CD with those of existing hybrid frameworks, which generally use one of two approaches. They couple separate codes through an external interface or embed a collisional $N$-body module directly into a hydrodynamic code.

An early example of the first approach is the Torch framework, which couples the grid-based adaptive mesh refinement MHD code Flash to collisional stellar dynamics through Amuse ([112, 111, 105]). Early versions used the fourth-order Hermite integrator Ph4 ([95, 137, 136]), together with the multiple-system module Multiples ([110]). Later versions replaced Ph4 with PeTar. Torch can simulate $\sim10^{6}\,M_{\odot}$ clouds and $\gtrsim10^{5}$ stars, and has been applied to SFE ([107]), runaway stars ([109]), primordial binaries ([19, 20]), and mass segregation ([109]). Another productive interface-based effort is the SIRIUS (SImulations Resolving IndividUal Stars) project, which developed the tree-direct hybrid SPH/N-body code Asura+Bridge (M. Fujii et al. 2007; M. S. Fujii et al. 2021a,b; M. Fujii et al. 2024). Applications include the formation history of the Orion Nebula Cluster ([40 b]), intermediate-mass black hole formation (M. Fujii et al. 2024), star cluster assembly ([76, 75]), and primordial binaries ([57]) during star cluster formation in molecular clouds. Compared with Arepo-CD, these frameworks can combine existing physics modules more flexibly, for example including modules for primordial binary stars. However, coupling through an external interface can add communication and synchronization overhead. Its effect on efficiency depends on the implementation and the problem. For example, Torch aggregates newly formed stars below a threshold mass to reduce the number of dynamical particles, sacrificing star-by-star sampling at the low-mass end. In contrast, Arepo-CD enables realistic star cluster formation simulations at moderate computational cost. One realization of the fiducial $10^{5}\,M_{\odot}$ GMC simulation requires only $\sim3000$ CPU hours,¹¹ including collisional dynamics and radiative transfer. In terms of scientific outcomes, Arepo-CD can produce similar results for SFE ([107]), runaway stars generated by dynamical ejection and subcluster mergers ([108]), and star cluster assembly ([76, 75]).

A complementary strategy is to embed collisional N-body methods directly within hydrodynamic codes. F. Dinnbier & S. Walch [33] and B. Zimmermann et al. [147] implemented a fourth-order Hermite integrator for sink particles in FLASH, although without full regularization of few-body systems. Another example is the use of the KETJU module [121, 92] within SPHGAL. N. Lahén et al. (2025b,a) extended KETJU to resolve collisional dynamics around massive stars with $m_\star> 8\,M_\odot$, applying it to idealized dwarf galaxies and to star cluster formation in a low-metallicity dwarf-galaxy merger. Such integrated implementations can reduce the overhead of coupling separate codes, but resolving close encounters only for a selected subset of stars provides an incomplete treatment of collisional dynamics. This approximation may be adequate for galaxy-scale simulations, but it is less suitable for the smaller-scale star cluster simulations considered here. In these systems, interactions among the full stellar population can be dynamically important.

More recent hydrodynamic codes with embedded collisional $N$-body solvers address this limitation by treating collisional dynamics for all stars. One example is Enzo-N (Y. Jo et al. 2024), which couples the collisional $N$-body code NBODY6++GPU to the cosmological MHD code Enzo. The framework was later rewritten and renamed Enzo-Abyss [17] (Y. Jo et al. 2026). This framework targets the formation of (nuclear) star clusters and their coevolution with host galaxies. It also targets high-redshift supermassive black hole seed formation through runaway stellar collisions in dense molecular clouds. Similarly, Y. Bernard et al. (2025) implemented a new collisional $N$-body scheme in Phantom. They simulated smaller molecular clouds than those considered here, with a higher gas mass resolution of $\Delta m_{\mathrm{gas}} \approx0.003\,M_\odot$. They also emphasized the importance of multiple realizations to reduce the impact of stochastic and chaotic variations in the $N$-body evolution. Enzo-Abyss and the Phantom implementation employ a two-level hierarchy combining a fourth-order integrator with a higher-order specialized integrator, whereas Arepo-CD uses a three-level hierarchy that combines second-order leapfrog, fourth-order FSI, and SDAR. This hierarchy applies different integration methods at different dynamical scales, providing greater flexibility for systems spanning a wide dynamical range. The computational cost of the present Arepo-CD simulation also makes multiple realizations and broader parameter surveys practical.

Future Extensions

As a proof-of-concept study, the present work focuses on the main numerical methods and their initial applications. However, the current implementation of Arepo-CD does not yet include all physical processes relevant to young star clusters. We discuss several possible extensions below.

An important future extension of the current framework is the treatment of stellar collisions and mergers. In very dense systems, repeated close encounters can lead to direct stellar collisions, whose remnants may evolve very differently from isolated stars [134, 135]. A fully self-consistent implementation would require a collision prescription and an evolutionary model for the merger products, including mixing and mass loss. This treatment is computationally demanding, although recent approaches based on machine learning may substantially reduce the cost [4, 48, 122]. A more practical first step is to enable the merger support already available in SDAR [140] and then couple it to an analytic prescription for subgrid merger outcomes [118]. This extension would allow Arepo-CD to explore the possibility of runaway stellar collisions and the self-consistent formation of very massive stars (VMSs). These simulations would start from a realistic GMC (M. S. Fujii et al. 2024) instead of specially constructed pure $N$-body ICs.

A further extension is the inclusion of primordial binaries. In the present simulations, compact binaries form only through dynamical interactions, demonstrating that few-body collisional evolution is active. However, real young stellar populations are born with substantial multiplicity. The initial binary fraction, mass-ratio distribution, orbital periods, and eccentricities can strongly affect early cluster heating, merger rates, and runaway production. A realistic primordial-binary module is also likely to increase the collision rate in dense embedded clusters and may promote the formation of VMSs [116]. An important next step is to implement this module consistently with both the IMF and the reservoir-based star formation scheme [30]. Encouragingly, pioneering work within the AREPO and RIGEL framework has already explored this direction (e.g. G. Liu et al. 2026 in prep.), using mature binary-pairing algorithms [80, 98, 19] that reproduce observed binary trends well [81].

A closely related extension is the inclusion of binary and multiple stellar evolution (BSE&MSE). The current Arepo-CD implementation follows only single-star evolution through the underlying RIGEL framework [27]. To model interacting binaries, mass transfer, common-envelope evolution, and compact-object binaries formed in dense young clusters, we will ultimately require coupling the subsystem integrator to a BSE or MSE package. This is technically challenging because the relevant dynamical and evolutionary timescales can become comparable in hierarchical few-body systems. Nevertheless, such a coupling will be essential for Arepo-CD to progress from early cluster dynamics to predictive studies of compact-object populations and merger channels. Established options include the classical BSE package [69] and more modern alternatives such as Sevn ([127]) and Binary_c ([71]). A longer-term option is the secular multiple-evolution code mse, which is explicitly designed to treat arbitrary hierarchical systems and their coupled dynamical and stellar evolution ([56]).

Beyond these planned extensions, the framework also points to several immediate scientific applications. The present GMC experiments highlight the roles of feedback timing, runaway production, and early binary formation in shaping embedded clusters. Extending these calculations to systematic surveys over cloud mass, density, metallicity, IMF sampling, and feedback prescriptions will be essential for direct comparison with observed young cluster populations and their scaling relations, such as the mass-radius relation of young star clusters ([14]). Over longer timescales, Arepo-CD should also enable studies of cluster survival, expansion, and tidal disruption in evolving galactic environments, and eventually in cosmological settings. In this sense, this study is a first step toward a predictive framework linking the small-scale physics of close stellar encounters to the large-scale processes of cluster formation and feedback in galaxies.

Summary and Conclusions

In this paper, we introduced Arepo-CD, a hybrid framework designed to follow star cluster formation in the regime where gas dynamics, stellar feedback, and collisional stellar dynamics operate simultaneously. Arepo-CD couples star-by-star collisional $N$-body dynamics to the moving-mesh MHD code Arepo. It uses unsoftened gravity for star-star interactions while retaining softened gravity for interactions involving gas and other collisionless components. Arepo-CD identifies compact stellar subsystems on-the-fly and evolves them with a hybrid integration scheme combining a fourth-order forward symplectic integrator with slow-down algorithmic regularization. The remaining system follows the native Arepo timestep hierarchy.

We first validated the method by comparing isolated star cluster simulations with those of the collisional $N$-body code PeTar. For the equal-mass model, Arepo-CD reproduces the evolution of the core and Lagrangian radii within the realization-to-realization scatter of PeTar. For clusters with realistic IMFs, the two codes also show close agreement in the radial density profiles and stellar velocity distributions. In particular, Arepo-CD reproduces the high-velocity escaper tail generated by strong collisional interactions, whereas simulations with softened stellar gravity remain more centrally concentrated and produce fewer fast escapers. These tests demonstrate that the hybrid integration scheme captures the main structural and kinematic signatures of collisional cluster evolution.

We then applied Arepo-CD to self-consistent GMC simulations coupled to the Rigel star-formation and feedback model. The collisional and collisionless runs exhibit very similar large-scale gas evolution, star formation rates, and integrated star formation efficiencies, indicating that turbulent collapse and stellar feedback primarily regulate global cloud evolution in the model studied here. Nevertheless, modest local differences develop in the stellar distribution and consequently in the spatial morphology of feedback. A much clearer difference appears in the stellar kinematics. Arepo-CD produces a pronounced high-velocity tail extending beyond 100 km/s and numerous runaway and walkaway stars, while the collisionless Rigel run produces no comparable runaway population. This demonstrates that strong few-body encounters can already become dynamically important during the embedded phase of cluster formation.

These results establish Arepo-CD as a practical and extensible framework that combines collisional stellar dynamics with hydrodynamic simulations of star cluster formation. Future extensions of the present implementation include stellar collisions and mergers, primordial binaries, and binary stellar evolution. Broader parameter surveys will enable systematic studies of runaway production, feedback redistribution, and the long-term evolution and survival of young star clusters in increasingly complex galactic environments.

Acknowledgments

We thank Volker Springel for giving us access to Arepo. CK is grateful to Chengzhe Li and Jiani Chu for useful discussions. HL is supported by the National Key R&D Program of China No. 2023YFB3002502, the National Natural Science Foundation of China under No. 12373006 and 12533004, and the China Manned Space Program with grant No. CMS-CST-2025-A10. LW acknowledges support from the National Natural Science Foundation of China through grants 12233013 and 12573041, the Fundamental Research Funds for the Central Universities, Sun Yat-sen University (2025QNPY04), the High-level Youth Talent Project (Provincial Financial Allocation) through the grant 2023YHSPY0706. YD is supported by the National Natural Science Foundation of China under No. 125B2057. EMH acknowledges support from the Tsinghua Shui Mu Scholarship and the Beijing Natural Science Foundation (IS25036).

Software: Arepo ([130]; [142]), MakeCloud ([50]), Matplotlib ([68]), McLuster ([84]; [139]), NumPy (S. van der Walt et al. 2011), PeTar ([138]), SDAR ([140]).