Issue 
A&A
Volume 649, May 2021



Article Number  A24  
Number of page(s)  11  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/202038784  
Published online  05 May 2021 
Introducing a new multiparticle collision method for the evolution of dense stellar systems
Crashtest Nbody simulations
^{1}
Dipartimento di Fisica e Astronomia & CSDC, Università di Firenze, Via G. Sansone 1, 50019 Sesto Fiorentino, Italy
email: pierfrancesco.dicintio@unifi.it
^{2}
INFN – Sezione di Firenze, Via G. Sansone 1, 50019 Sesto Fiorentino, Italy
^{3}
INAF, Osservatorio Astronomico di Padova, Vicolo dell’Osservatorio 5, 35122 Padova, Italy
^{4}
Center for Astro, Particle and Planetary Physics (CAP 3), New York University, New York, Abu Dhabi, UAE
email: mp5757@nyu.edu
^{5}
INFN – Sezione di Padova, Via Marzolo 8, 35131 Padova, Italy
^{6}
Department of Astronomy & Center for Galaxy Evolution Research, Yonsei University, Seoul 03722, Republic of Korea
email: sjyoon0691@yonsei.ac.kr
Received:
29
June
2020
Accepted:
18
November
2020
Context. Stellar systems are broadly divided into collisional and noncollisional categories. While the latter are largeN systems with long relaxation timescales and can be simulated disregarding twobody interactions, either computationally expensive direct Nbody simulations or approximate schemes are required to properly model the former. Large globular clusters and nuclear star clusters, with relaxation timescales of the order of a Hubble time, are small enough to display some collisional behaviour and big enough to be impossible to simulate with direct Nbody codes and current hardware.
Aims. We aim to introduce a new method to simulate collisional stellar systems and validate it by comparison with direct Nbody codes on smallN simulations.
Methods. The MultiParticle Collision for Dense Stellar Systems (MPCDSS) code is a new code for evolving stellar systems with the multiparticle collision method. Such a method amounts to a stochastic collision rule that makes it possible to conserve the exact energy and momentum over a cluster of particles experiencing the collision. The code complexity scales with N log N in the number of particles. Unlike Monte Carlo codes, MPCDSS can easily model asymmetric, nonhomogeneous, unrelaxed, and rotating systems, while allowing us to follow the orbits of individual stars.
Results. We evolved small (N = 3.2 × 10^{4}) star clusters with MPCDSS and with the directsummation code NBODY6, finding a similar evolution of key indicators. We then simulated different initial conditions in the 10^{4} − 10^{6} star range.
Conclusions. MPCDSS bridges the gap between small collisional systems that can be simulated with direct Nbody codes and large noncollisional systems. In principle, MPCDSS allows us to simulate globular clusters such as Ω Centauri and M 54, and even nuclear star clusters, which is beyond the limits of current direct Nbody codes in terms of the number of particles.
Key words: methods: numerical / Galaxy: bulge / globular clusters: general / galaxies: dwarf
© ESO 2021
1. Introduction
Our understanding of the formation and dynamical evolution of dense stellar systems, such as globular clusters (hereafter GCs) and nuclear star clusters (hereafter NSCs), has a crucial impact on galactic archaeology (see e.g., Belokurov et al. 2018; Myeong et al. 2018, 2019; Massari et al. 2019; Ibata et al. 2019; Di Matteo et al. 2019; Chung et al. 2019; Grand et al. 2020), multimessenger astronomy (where it allows us to better constrain compact object mergers; see e.g., Belczynski et al. 2002; Banerjee et al. 2010; Bae et al. 2014; Ziosi et al. 2014; Breivik et al. 2016; Rodriguez et al. 2016; Hurley et al. 2016; Askar et al. 2017; Chatterjee et al. 2017; Arca Sedda et al. 2018; Kremer et al. 2018; Di Carlo et al. 2019; Bouffanais et al. 2019; Rastello et al. 2019; Antonini & Gieles 2020), cosmology and supermassive black hole science (with star clusters acting both as nurseries of intermediatemass black hole seeds and a delivery mechanism to the galactic centres; see e.g., CapuzzoDolcetta 1993; Ebisuzaki et al. 2001; Portegies Zwart et al. 2004, 2006; CapuzzoDolcetta & Miocchi 2008; Antonini et al. 2012; MastrobuonoBattisti et al. 2014; ArcaSedda & CapuzzoDolcetta 2014; Askar et al. 2021), and even stellar astrophysics (as clusters are key to the formation of stellar exotica; see e.g., Fabian et al. 1975; Bailyn 1995; Portegies Zwart et al. 2001, 2010; Fregeau et al. 2004; Verbunt & Lewin 2006; Leigh et al. 2007; Pasquato et al. 2014; van den Berg 2019; Wang et al. 2020).
However, modelling selfgravitating Nbody systems with a realistic number of stellar particles (in some cases well above 10^{6}) over several relaxation times is extremely challenging in terms of computational resources. This is due to the superquadratic scaling of complexity with the number of particles in direct summation codes (Makino & Hut 1988; Aarseth 2003). Current stateoftheart direct Nbody simulations with 10^{6} particles need several months of computer time – even on dedicated accelerator graphic processing unit (GPU) clusters – to follow the evolution of typical globular clusters (Wang et al. 2016). This is an issue especially when simulating collisional systems, where the effects of relaxation driven by twobody interactions cannot be neglected.
In Fig. 1, we show that a large fraction of the Milky Way globular clusters is both in the collisional regime and contains a sufficiently high number of stars to make detailed modelling based on direct Nbody simulations infeasible, especially when the need to obtain a significant number of realisations of the same system is taken into account.
Fig. 1. Globular star clusters of the Milky Way and satellites from McLaughlin & van der Marel (2005). Log halfmass relaxation time (in years) is plotted against log total mass (in Solar mass). Globular clusters Ω Centauri and M 54 are shown as a magenta square and cyan triangle, respectively. Assuming a mean stellar mass of 0.5 M_{⊙}, star clusters with a mass above 10^{5} M_{⊙} (shaded in light blue) contain 2 × 10^{5} stars and can be simulated only with great computational effort, and they are limited to a handful of realisations. This is especially true if a realistic binary fraction is included, at variance with smaller systems (grey shaded area) that are well within the capabilities of current direct Nbody codes. Star clusters with a relaxation time over the typical globular cluster age (≈10 Gyr) can be regarded as collisionless and are shaded in peach and dark green. According to this definition, Ω Centauri and M 54 are only slightly collisional, but clearly not accessible to modelling through direct Nbody simulations. 
Several approximated alternatives to the direct Nbody approach^{1} that do not share its prohibitive computational cost exist (see Heggie 2016 for an excellent review). The family of socalled Monte Carlo codes, which essentially solve the FokkerPlanck equation, is perhaps the most successful among these (Hénon 1971a,b, 1975; Stodolkiewicz 1982, 1986; Joshi et al. 2000; Freitag & Benz 2001, 2002; Giersz 2001, 2006; Pattabiraman et al. 2013; Giersz et al. 2013; Hypki & Giersz 2013; Pijloo et al. 2015; Rodriguez et al. 2018; Sollima & Ferraro 2019). Monte Carlo simulations, however, are generally limited to spherically symmetric systems, with the notable exception of the code developed by Vasiliev (2015), which is unfortunately not in widespread use. Among other issues, this limitation was shown to lead to discrepancies between direct Nbody and Monte Carlo (or any FokkerPlanck solver that assumes spherical symmetry) in the presence of an external tidal field such as the Galactic one (Takahashi & Portegies Zwart 2000), with the notable exception of the scheme by Sollima & Mastrobuono Battisti (2014), which shows remarkably good agreement with direct Nbody simulations with N ≈ 2 × 10^{4} when applied to GCs orbiting in a fixed pointlike galactic potential.
In this first paper of a series we introduce a new simulation scheme, the MultiParticle Collision for Dense Stellar Systems (hereafter MPCDSS) code, which combines an essentially linear scaling of computational complexity in the number of particles with the ability to model configurations with arbitrary geometries. Here, we introduce the structure of the code and present an initial series of tests on the dynamical evolution of GCs, without including the stellar evolution modules, paving the way to the application to the two most massive star clusters in the Milky Way^{2}: M 54 and Ω Centauri. Both clusters show a spread in metallicity (Sarajedini & Layden 1995; Lee et al. 1999), hinting at a nontrivial dynamical history that possibly includes one or more mergers (see e.g., Cole et al. 2017; AlfaroCuello et al. 2020a,b) that may still affect presentday observable properties (AmaroSeoane et al. 2013; Pasquato & Chung 2016). The clusters have large masses, thus going beyond the limit of what is currently modellable with ‘honest’ direct Nbody simulations.
This paper is structured as follows: in Sect. 2, we introduce the numerical methods used in MPCDSS to compute the gravitational field, treat the collisions and propagate the simulations particle trajectories, and discuss the efficiency of our implementation. In Sect. 3, we compare a set of test simulations of collisional evaporation of smaller GCs using MPCDSS and the stateoftheart direct Nbody code NBODY6 (Aarseth 2003; Nitadori & Aarseth 2012). In Sect. 4, we present the results of numerical simulations of core collapse and mass segregation. In Sect. 5, we discuss our findings, and, we summarize our results.
2. Overview of the numerical method
2.1. Stochastic collisions: The multiparticle collision method
In our numerical code, we resolve the collisional interactions between stars using the socalled multiparticle collision method (hereafter MPC). Originally, MPC was introduced by Malevanets & Kapral (1999, 2004) in the context of numerical hydrodynamics for the simulation of mesoscopic fluids (e.g., polymers in solution, colloidal fluids). It has been shown that the method yields Galileaninvariant dynamics, that the NavierStokes equations are recovered in the continuum limit, and that relaxation towards thermodynamical equilibrium is correctly modelled (see Gompper et al. 2009 for a detailed review).
Recently, the MPC techniques have also been used in plasma physics to address heat transport problems in reduced models in 1D (Di Cintio et al. 2015; Ciraolo et al. 2018; Lepri et al. 2019) and 2D (Di Cintio et al. 2017a,b, 2018).
The MPC scheme alternates a streaming step (corresponding to noncollisional evolution) and a collision step. In 3D space the collision step amounts to a rotation of the particle’s velocity vectors in the centre of mass frame of each cell^{3} onto which the simulation domain has been partitioned.
At the beginning of the collision step, the code evaluates the centreofmass (c.o.m.) velocity in every cell:
and the relative velocities δv_{j} = v_{j} − u_{i}. Then, for each cell a random axis R_{i} and rotation angle α_{i} are sampled from uniform distributions. At this point, the vectors δv_{j} are rotated around R_{i} of α_{i} and then converted back to the simulation frame, so that for the jth particle in cell i the new velocity reads
where δv_{j, ⊥} and δv_{j, ∥} are the relative velocity components perpendicular and parallel to R_{i}, respectively.
Such an operation conserves exactly the total kinetic energy K_{i} and the three components of the momentum P_{i} in cell i (see e.g., Ryder 2005; Di Cintio et al. 2017a for the rigorous proof). By introducing a constraint on the rotation angles α_{i}, we conserve a component of the angular momentum vector of the cell L_{i} by defining α_{i} such that
where
In the formulae above, r_{j} are the particles’ position vectors, and the notation [x]_{z} means that one is taking (without loss of generality) the component of the vector x parallel to the z axis of the simulation’s coordinate system, so that the z component of the cell angular momentum is conserved.
We note that, for strictly 2D systems, Eq. (2) becomes
where 𝒢_{αi, i} is now the 2D rotation matrix of an angle α_{i} that, if chosen according to Eqs. (3) and (4), ensures the conservation of the scalar angular momentum (see Di Cintio et al. 2017a) in addition to K_{i} and P_{i}. We also note that the conservation of the total angular momentum can be achieved even in 3D systems by choosing R_{i} to be parallel to the direction of the cell’s angular momentum vector L_{i} and taking with the notation above. For the simulations presented here, we limited ourselves to the standard rotation scheme with only one component of the total angular momentum conserved, as it is much less timeconsuming not having to determine the direction of the angular momentum (pseudo)vector cell by cell.
As in GCs the collision frequency strongly depends on the local values of the stellar density and velocity dispersion, we conditioned the MPC step to a celldependent probability accounting for the local degree of collisionality. We first defined the celldependent MPC probability as
where Δt is the timestep, the mean stellar number density, and σ_{i} the average particle mass and the velocity dispersion in the cell, respectively, β is a dimensionless constant fixed to 2N_{c}, and Erf(x) the standard error function. The celldependent Coulomb logarithm is defined as
with r_{s} being the typical scale length of the system.
Once Eq. (6) is evaluated in each cell, a random number p_{*i} is sampled from a uniform distribution in the interval [0,1] and the multiparticle collision is applied for all cells for which p_{*i} ≤ p_{i}. By doing so, particles in cells with higher collision frequency are more likely to be prone to a MPC step.
We note that applying a MPC scheme puts an additional constraint on the choice of the simulation time step Δt. In particular, if the latter is too short compared to a given cellcrossing time (i.e. the time it takes to a test particle moving at the local mean speed to cross the cell), the collisions could be overestimated. Viceversa, if Δt is too large, some fast particles might not spend a sufficient amount of time in cells at large collisionality (i.e. with local large density), virtually never experiencing collisions.
2.2. Deterministic dynamics
We simulated the collective dynamics of the systems using a standard particlemesh scheme (see e.g., Hockney & Eastwood 1988) that solves the Poisson equation:
on a spherical grid in polar coordinates N_{r} × N_{ϑ} × N_{φ} and interpolates ∇Φ at each particle position r_{i}.
The equations of motion between two MPC steps were solved with a standard second order leapfrog scheme with a fixed timestep (see e.g., McLachlan & Atela 1992) of the order of 10^{−2} in units of the system’s initial crossing time t_{dyn} (see Eq. (12) below). In the preliminary simulations presented in this work, in order to further speed up the calculations, instead of solving Eq. (8), we only evaluated the radial component of the gravitational field so that, in practice, the equations of motion became
where M(r_{i}) is the mass within the particle’s radial coordinate r_{i}. By doing so, when needed, the potential Φ(r_{i}) on particle i can be obtained as
after having sorted the radial coordinates of all particles (see e.g., Pattabiraman et al. 2013; Rodriguez et al. 2018) as in standard Monte Carlo codes. With such assumptions, the initial spherical symmetry of the model would be preserved, as no radial orbit instability can take place (Ciotti et al. 2007), and, if the collision step were deactivated, each particle orbit would retain its original plane, and all angular momentum vectors J_{i} = m_{i}r_{i} × v_{i} would be conserved individually.
3. Comparison with direct Nbody
Introducing a somewhat stochastic mechanism of scattering between individual particle orbits via the MPC scheme might at first seem to be inducing dramatic differences in the dynamical behaviour of orbit evolved with MPCDSS with respect to standard direct Nbody integrators.
3.1. Conservation laws and orbital structure
As a first test of reliability of MPC simulations, we evolved small systems of 10 000 ≤ N ≤ 32 000 equal mass particles with the direct Nbody code NBODY6 and with MPCDSS and compared the evolution of the virial ratio −2K/U (where K is the total kinetic energy and U the total potential energy), and the total energy and angular momentum fluctuations δE = (E − E_{0}) and δL = (L − L_{0}) (where E_{0} and L_{0} are the initial values of the energy and norm of the angular momentum, respectively), as resolved by the two schemes for identical isotropic equilibrium initial conditions. In addition, we also studied the structure of the individual tracer particle orbits under both numerical methods.
The models used for the two simulations sets were set up as follows. We considered an isolated spherical isotropic model with the following Plummer (1911) density distribution:
with total mass M and scale radius r_{s}, related to the halfmass radius by r_{half} ≈ 1.3 r_{s}. In both the Nbody and MPC cases, the systems were evolved up to 1000 dynamical times, defined by
with a fixed timestep Δt = t_{dyn}/100 and neglected stellar evolution and formation of binaries (such that each simulation particle represents an individual star retaining its mass throughout the whole run). With such a choice of Δt, we retain a good energy conservation while avoiding the overestimation of the stochastic collisions in a given cell.
As an example of the conservation of an equilibrium state, in Fig. 2 we show the evolution of the virial ratio as well as the total energy and angular momentum fluctuations, δE and δL, for an equilibrium isotropic Plummer (1911) model with 2 × 10^{4} particles evolved with the direct Nbody code and the two versions of MPCDSS with and without the angular momentum conserving scheme. We only carried out the simulation up to 500 t_{dyn} so that the collisions would not affect the model’s properties much. In all cases, we used the same fixed timestep Δt and a secondorder propagation scheme.
Fig. 2. Evolution of the the virial ratio −2K/U (top panel), the fluctuations of the total energy δE in units of the initial energy E_{0} (middle panel), and the norm of the total angular momentum δL (bottom panel) for the same isotropic Plummer initial conditions with N = 2 × 10^{4} particles, evolved with the direct Nbody code (black lines), the exact angular momentum preserving MPC scheme (light blue lines), and the faster MPC rotation with a random axis (purple lines). The initial conditions have a slight angular momentum due to the randomised initialisation procedure for stellar velocities. 
Remarkably, in the two MPC simulations, the virial equilibrium is preserved with smaller fluctuations with respect to the direct Nbody integration. This is due to the fact that the enforced spherical symmetry and the smoother gridbased potential in the MPCDSS simulations generally induce smaller force fluctuations on particles (for a given choice of N) with respect to a direct force calculation. Accordingly, the total energy only shows fluctuations δE/E_{0} as small as 10^{−3} for the MPC simulations at variance with the direct Nbody simulations where δE is already of the order of 10^{−2} at early times. When using the MPC scheme with enforced angular momentum conservation, for this particular choice of initial conditions and simulation setup, its fluctuations δL/L_{0} are of the order of 10^{−4}, while they are of the order of 10^{−3} for the direct Nbody simulation (see the inset in the lower panel of Fig. 2). When the local conservation of angular momentum is not implemented in the MPC step, the total angular momentum might experience wild fluctuations δL/L_{0} of the order of 0.5 the smaller the initial value of the angular momentum L_{0}. Remarkably, breaking the conservation of angular momentum does not affect the energy conservation nor the virial equilibrium.
Concering the behaviour of individual particle orbits, using the MPC step to resolve collisional processes surprisingly does not qualitatively alter the behaviour of the orbits themselves with respect to standard Nbody simulations. In the upper panels of Fig. 3, we show the projections on the x − y plane of two orbits in a Plummer model with N = 2 × 10^{4} integrated with a direct Nbody code (panels a and b, black lines) and with MPCDSS (panels c and d, red lines). Always confined within less than two halfmass radii, the particles experience several ‘close encounters’ in both cases, thus subjecting them to dramatic changes in orbital inclination and precession frequency. Direct Nbody and MPC dynamics results in a large degree of phasespace ‘diffusion’ with particles exploring the whole energetically accessible region as shown in the lower panels of Fig. 3.
Fig. 3. Projections on the x − y plane of two tracer orbits starting from the same initial conditions propagated in a Plummer system with N = 2 × 10^{4} equal mass particles, evolved with an Nbody code (panels a and b, black lines) and MPCDSS (panels c and d, red lines). Bottom panels (e–h): corresponding phasespace sections in the r − v_{r} subspace. 
Moreover, following Di Cintio & Casetti (2019), we also studied the Fourier spectra of the radial coordinate r for individual particle orbits for a broad range of initial conditions. In general, the stochastic collision rule does not significantly alter the structure of the spectrum of a given orbit obtained from the same initial condition, with respect to a direct Nbody evolution. In Fig. 4, we show the modulus squared r^{*}(ω)^{2} of the radial coordinate r for a particle propagated in the same Plummer model with N = 2 × 10^{4} with the two simulation approaches, corresponding to panels a and c of Fig. 3. Remarkably, the structure of the fundamental frequency (and a large fraction of higher harmonics at larger values of ω) is preserved, thus leading to speculation that MPC dynamics can be sufficiently trusted even for larger systems, as the introduction of a stochastisation of particle velocities does affect the collective behaviour of the models (even for a rather small N such as 2 × 10^{4}). In addition, for the same orbits of Fig. 4, in Fig. 5 we show the probability density function of their x coordinate 𝒫(x) (i.e. the distribution of positions along x attained by the tracer particles within a given time interval) up to t = 1000 t_{dyn}. The two distributions match remarkably well, thus confirming the overall agreement between the individual orbital structure in the two simulation methods. We note that a similar comparison to Nbody orbits has also been made for massive tracer particle orbits propagated semianalitically with FokkerPlanck schemes (see e.g., Chatterjee et al. 2002a,b, 2003), or with direct integration of the Langevin equation (Di Cintio et al. 2020) in N = 10^{5} equal mass Plummer models. In general, at least for the models considered in these studies, the distribution of the coordinates 𝒫(x) evaluated from the numerical solution of the FokkerPlanck equation has the best agreement with the results of the direct Nbody simulations, with the Langevin approach being strongly dependent on the form of the noise term. Our MPC approach places itself in between the two methods.
Fig. 4. Fourier spectra of the radial coordinate r for typical orbits extracted from an Nbody simulation (red curve) and an MPC simulation (black curve) for a system of N = 2 × 10^{4} equal mass particles. 
Fig. 5. Probability density function of the x coordinate 𝒫(x) for a tracer particle starting with the same initial condition in a N = 2 × 10^{4} Plummer model evolved with the direct Nbody (black circles) and the MPC (red squares) methods. 
We note that using the MPC scheme results in a dramatic reduction of the computational cost of the simulation while retaining sufficiently reliable results. For example, the evolution of an N = 32 000 cluster up to 2000 t_{dyn} takes a few days on a dedicated GPU workstation for the direct Nbody simulation without stellar evolution and adaptive time step. When using an N_{c} = 32 × 16 × 16 grid in polar coordinates to perform the collisions and an N_{r} = 2000 radial grid to evaluate the gravitational field in monopole approximation (see Eq. (9)), our MPC simulations with fixed Δt take roughly one hour on a single core of a 64 bit INTEL® machine.
3.2. Core collapse with mass spectrum
The effect of multiple mass populations on the dynamical evolution of GCs is of prime importance. The present implementation of our MPC code offers the interesting chance to study multimass systems without adding extra computational complexity. We performed a set of additional tests with NBODY6 and MPCDSS simulating the evolution up to and after core collapse of Plummer models with a mass spectrum.
For the sake of simplicity, instead of using the multislope Kroupa (2002) mass function, in this work the particle masses m_{i} were extracted from a pure powerlaw mass function of the form
where α > 0, and the normalisation constant C depends on the minimumtomaximum mass ratio ℛ = m_{min}/m_{max} so that .
We ran simulations with N = 32 000 for a range of α spanning from 0.6 to 3.0 by intervals of 0.1, and in the case of the direct Nbody simulations we ran 10 different realisations (with a different seed for the initial conditions) for each value of α. In all cases with mass spectra discussed in this work, we fix ℛ to 10^{−3}.
In both the Nbody and MPC cases, the systems were evolved for 10^{4} dynamical times, corresponding to roughly 20 twobody (collisional) relaxation times of a model with the same total mass and number of equal mass particles, given by
Again, as a rule, in all sets of MPCDSS simulations we fixed the timestep Δt to t_{dyn}/100. We neglected stellar evolution and the contribution of binaries^{4}. For the two simulation sets (Nbody and MPC), we extracted and compared the number of escapers N_{esc} (defined as the number of particles at r > 17 r_{s} with positive energy; see e.g., Fukushige & Heggie 2000), the central density and velocity dispersion ρ_{0} and σ_{0}, as well as the mass function in the core as a function of time.
In Fig. 6, we present the time evolution of the number of escapers N_{esc} for choices of α in Eq. (13). We find that in all cases the MPC evolution (solid line) recovers the quasilinear trend of N_{esc} with time. However, some discrepancies between MPC and Nbody simulations are observed, in particular at low values of α (i.e. ‘flatter’ mass spectrum). We interpret such a difference as the effect of strong twobody collisions between light and heavy particles, resolved in a direct Nbody code but somewhat smearedout in a multiparticle collision code. In practice, in models with a mass spectrum with a larger fraction of heavy particles, it is more likely for the lighter ones to be ejected with positive energy when experiencing close encounters with heavier stars as compared with systems with mass spectra strongly peaked at low masses. Since the MPC step (even in the presence of a mass spectrum) simulates multiple interparticle collisions with a single move acting on all particles of the cell, the contribution of few but strong lightheavy particle encounters is obviously underestimated.
Fig. 6. Number of escapers N_{esc} as a function of time in units of the dynamical time t_{dyn} for a Plummer model with N = 32 000 and powerlaw mass spectrum with, from left to right, α = 1.5, 2.0, 2.3, 2.5, and 3.0 when evolved with NBODY6 (dashed lines) and with MPCDSS (solid lines). 
For the best agreement case, represented by the simulations with α = 2.3 (and corresponding to a Salpeter 1955 mass function), in Fig. 7 we show the evolution of the central density ρ_{0} and central velocity dispersion σ_{0} defined within the Lagrangian radius enclosing 8% of the total mass M. In this case, the evolution of both quantities with the MPC code (squares) matches that obtained using the NBODY6 (circles) remarkably well.
Fig. 7. Evolution of the central density ρ_{0} (upper panel) and central velocity dispersion σ_{0} (lower panel) for the Nbody simulation (empty symbols) and the MPC simulation (filled symbols) for the model with α = 2.3 in Fig. 6. The symbol sizes are of the order of the mean (Poissonian) error of particle counts in both plots. 
In addition, for each simulation of the two sets we took the time of core collapse t_{cc} as the time at which the minimum value reached by r_{2%} (i.e. the 3D Lagrangian radius enclosing 2% of the total mass). Figure 8 shows the time of core collapse as a function of the initial massfunction powerlaw exponent α. As expected, simulations starting with a shallower mass function have heavier particles, slowing down the core collapse. Thus, lowα runs take longer to reach corecollapse. Interestingly, for large values of α (say ≳2.6), t_{cc} starts increasing again in both MPC and direct Nbody simulations. We speculate that such curious behaviour might be due to the less efficient dynamical friction in models dominated by lowmass particles (see Ciotti 2010) resulting in larger sinking timescales for particles sitting in the largemass part of the spectrum. To guide the eye, we fitted the relation between α and t_{cc} using a secondorder polynomial (solid green line). We performed additional simulations (to be published elsewhere) with different values of N, α, and ℛ, and the nonmonotonic trend of t_{cc} with α appeared to be robust. Interestingly, the time of core collapse estimated with MPC simulations appears to be systematically larger than its counterpart for direct Nbody simulations for values of α larger than 2.3. We speculate that the reason for this discrepancy lies in the absence of a binary formation mechanism in the implementation of MPCDSSE used here. This, combined with the less efficient dynamical friction results in the time at which the collapse inverts being delayed, thus yielding a larger t_{cc}.
Fig. 8. Time of core collapse t_{cc} in units of the dynamical time t_{dyn} as a function of the mass spectrum exponent α for MPC (empty red squares) and Nbody (small black circles) simulations. A parabolic fit to the Nbody simulations (thick green line) is superimposed. The MPC result is generally within the range spanned by the direct Nbody realisations, except for a few values of α at the high end, where core collapse happens earlier in direct Nbody simulations. 
4. Further numerical experiments and results
We performed a broad spectrum of numerical experiments to determine the range of applicability of our newly introduced simulation method. First of all, we investigated the process of core collapse of singlecomponent models starting with isotropic Plummer density profiles in a broad range of system sizes spanning from 10^{3} to 10^{6}. For this set of numerical experiments, we followed the evolution of the 3D and associated Lagrangian radii containing different fractions of the total number of particles (or mass) between 2% and 90%.
In line with expectations, we find from the MPC simulations that for large values of N, the time of core collapse becomes asymptotically larger (in units of the dynamical time t_{dyn}), as shown in the left panel of Fig. 9 for the equal masses case, with N ranging from 10^{4} to 10^{6}. Remarkably, for low values of N, there seems to be a somewhat nonmonotonic trend of t_{cc} with N. We performed additional MPC simulations with N as small as 10^{4} using different realisations of the initial condition with choices of timestep and grid size, and such a trend persists.
Fig. 9. Relations of the number of particles N versus the time of core collapse t_{cc} in units of the dynamical time (left panel) and relaxation time (right panel) for singlecomponent Plummer models. 
When expressing t_{cc} in units of the collisional relaxation time t_{2b} (cf. Eq. (14)), the picture is inverted and large N systems reach core collapse earlier in units of their intrinsic t_{2b} (see right panel of the same figure). As N increases to larger and larger values, the systems start to approach the socalled collisionless limit, where the effects of collisions become more and more negligible, and therefore the process of core collapse is driven mainly by collective instabilities that take place on increasingly large timescales. On the other hand, as N increases, t_{2b} also increases as N/log N, resulting in a decreasing trend of t_{cc}/t_{2b} with N.
In Fig. 10, we show the radial density profile at different times between 1 t_{2b} and 18 t_{2b} for the N = 10^{5} case. It is evident, as at around 2 t_{2b} (corresponding roughly to 2340 t_{dyn} for this value of N) the density has already departed significantly from the initial Plummer profile (marked in figure by the thin black line). Remarkably, at later times, the inner part of the density slope approaches the r^{−2.23} trend (heavy dashed line) as predicted by Cohn (1980, see also Heggie & Stevenson 1988), and it is in nice agreement with the Monte Carlo simulations by Joshi et al. (2000) and Pattabiraman et al. (2013). For times larger than roughly 8 t_{2b}, the reexpansion of the outer regions becomes evident, as can also be seen from the evolution of the Lagrangian radii in Fig. 11. Surprisingly, even if in the current implementation of MPCDSS there is no explicit treatment of the binaries, the evolution of the density profile towards and beyond the core collapse matches remarkably well with its counterpart obtained with different numerical methods and incorporating a systematic treatment of binaries. In particular, the fact that the critical r^{−2.23} density slope is recovered in MPC simulations (even for systems as small as N = 10^{4} and as large as N = 10^{6}, not shown here) suggests that it is not strongly related to the ‘thermostat effect’ of binary formation, but rather to the evaporation of particles from the contracting core.
Fig. 10. 3D density profile for a model with N = 10^{5} equal mass particles and initial Plummer density distribution (thin black solid line) at t = 1, 2, 3, 6, 8, 10, and 18 twobody relaxation times t_{2b} (coloured lines). The heavy dashed line marks the theoretical r^{−2.23} profile. 
We find a corecollapse time t_{cc} (i.e. the time at which the central part of the cluster reaches the highest concentration that we measure here as the minimum attained by the Lagrangian radius containing the 2% of the simulation particles) of about 10 t_{2b} (indicated by the vertical line in Fig. 11), which is in rather good agreement with the Monte Carlo simulations of Joshi et al. (2000) and Hurley & Shara (2012), and the Nbody simulations of Küpper et al. (2008), finding values between 10 and 15 t_{2b} for models with initial conditions analogous to the ones used in our simulations (i.e. Plummer profile, N = 10^{5}, and no mass spectrum).
Fig. 11. Evolution of the 3D Lagrangian radii enclosing 2%, 5%, 10%, 50%, and 90% of the total number of particles N for the same model in Fig. 10. The vertical dashed line marks the system’s corecollapse time t_{cc} ≈ 10 t_{2b} ≈ 12 000 t_{dyn}. 
Moreover, we also performed additional simulations for Plummer models with different values of N and mass spectra, finding that, surprisingly, for the models with α in the range between 1.5 and 3, the asymptotic slope of the density profile in the inner regions matches better to the predicted r^{−2.23} trend, as shown in Fig. 12 for the α = 1.5 and α = 2 cases with N = 10^{5}.
Fig. 12. Evolution of the 3D density profile for two models with N = 10^{5} and initial Plummer density distribution (thin solid line) and mass spectra with α = 1.5 (left panel) and 2.0 (right panel). As in Fig. 10, the heavy dashed line marks the theoretical r^{−2.23} profile. 
We observe that, in general, for a fixed number of particles and a total mass M, models with a mass spectra reach the core collapse faster in units of t_{dyn} than the associated equal mass case. This can be seen from Fig. 13, where we show the Lagrangian radii enclosing the same fractions of the total number of particles as in Fig. 11, but for two models with mass spectrum with exponents α = 1.5 and 2. In both cases, the core collapse time t_{cc} is well below 10^{3} t_{dyn}, being roughly 2000 t_{dyn} for α = 1.5 and 773 t_{dyn} for α = 2.0 as marked by the vertical lines in the figure.
Fig. 13. Evolution of the 3D Lagrangian radii enclosing 2%, 5%, 10%, 50%, and 90% of the total number of particles N = 10^{5} for the two models shown in Fig. 12. The solid and dotteddashed lines refer the α = 1.5 and 2 cases, respectively. The two vertical dashed lines mark the corecollapse times for the models: t_{cc} ≈ 2000 t_{dyn} for α = 1.5 and t_{cc} ≈ 773 t_{dyn} for α = 2.0. 
Since the timedependent radii enclosing a given fraction of the total mass M or number of particles N are not the same quantity for a model with different species, we evaluated both types of ‘Lagrangian radii’ for some of the models with mass spectra. As expected, since as result of the more efficient dynamical friction, more massive stars tend to accumulate to the centre of the system, the Lagrangian radii computed for a given percentage of the mass of the model attain systematically smaller values than those evaluated for the same percentage of the total number of particles instead. However, the estimated corecollapse times do not differ significantly depending on the two choices, independently of the number of particles in the simulation.
5. Discussion, conclusions, and future prospects
We introduce a new code for simulating the collisional evolution of dense stellar systems using the MPC approach. Our code is characterised by an N log N scaling in the number of stars, which makes it suitable for simulating massive GCs and the Milky Way NSC. These systems are, for the foreseeable future, well beyond the reach of directsummation Nbody codes, because the latter scale quadratically with the number of stars. The MPC method is based on alternating streaming steps (where stars evolve in the smooth gravitational potential of the whole star system) and collision steps, which are intended to model the relaxation effects induced by stellar encounters. Collision steps are cellbased rotations of the stellar velocity vectors that, by construction, conserve mass, momentum, and energy^{5}. The MPC approach avoids the intricacies of two and multiplebody encounters that require indirect Nbody code techniques such as softening or KustaanheimoStiefel regularisation (see e.g., Mikkola 2008) in order to treat the singularity of the 1/r gravitational potential^{6} while retaining the relaxation effects of the encounters. Eliminating the need to compute all pairwise forces, the MPC approach results in much lower algorithmic complexity with the number of stars, without losing the ability to correctly recover the longterm evolution of stellar systems.
Compared to Monte Carlo approaches, our code has the advantage of easily simulating any geometry as the MPC method has no apparent dependence on the shape of the individual cells or the overall grid structure, while most Monte Carlo codes are confined to highly symmetric configurations. We can thus simulate rotating^{7}, merging, or tidally disrupted star clusters with no additional effort and no loss in accuracy with respect to spherically symmetric systems.
In this paper, we present our MPC code and a number of test simulations, showing that the total energy and angular momentum of an isolated simulation are conserved. We also calculated the virial ratio (kinetic over potential energy), which is also conserved with remarkable accuracy. Finally, we validated our code by comparison to direct Nbody simulations of star clusters. We find that the evolution of the central density, central velocity dispersion, and number of escapers as a function of time in MPC simulations closely follows that in direct Nbody simulations over a wide a range of stellar mass spectra. Additionally, the time at which core collapse is reached is also in good agreement with both direct Nbody and theoretical analytical calculations.
In the future, we plan to add an interface with a community stellar evolution module such as single stellar evolution (SSE; Hurley et al. 2000), or the more recent SEVN (Spera & Mapelli 2017), and to introduce one or more schemes to simulate binary stars, such as tracked particles with internal dynamics in the same spirit of plasma codes including ionisation and recombination (see e.g., Di Cintio et al. 2013), or a nested direct integrator for a subset of particles (e.g., Fregeau 2012) coupled with a binary evolution code (Hurley et al. 2002; Giacobbo et al. 2018). The aim will be to study the dynamics of compact objects (e.g., Mapelli 2016; Rastello et al. 2020) and other phenomena hinging critically on binary modeling, such as blue stragglers (Miocchi et al. 2015; Pasquato et al. 2018; Pasquato & Di Cintio 2020). We will then run a large set of simulations that we plan to release publicly along with the fully parallelised version of the code, including simulations intended to model specific objects, such as Ω Centauri and M 54. Finally, we will address the black hole retention problem in star clusters, focusing in particular on the fate of intermediatemass black holes in the Galactic NSC.
We note that direct Nbody solutions are themselves not exact, due to the chaotic nature of the problem (see e.g., Di Cintio & Casetti 2019, 2020 and references therein) and the finite precision of the numerics involved (see e.g., Breen et al. 2020, which also proposes a creative alternative simulation scheme).
With the exception of the NSC (Walcher et al. 2005; Misgeld & Hilker 2011; Neumayer 2017), which is in principle also amenable to simulation with MPCDSS and has also been studied with direct Nbody simulations (Agarwal & Milosavljević 2011; Perets & MastrobuonoBattisti 2014) in the context of the socalled repeated accretion scenario (Antonini et al. 2012; ArcaSedda & CapuzzoDolcetta 2014).
Acknowledgments
This project has received funding from the European Union’s Horizon 2020 research and innovation program under the Marie SkłodowskaCurie grant agreement No. 664931. This material is based upon work supported by Tamkeen under the NYU Abu Dhabi Research Institute grant CAP3. P.F.D.C. wishes to thank the financing from MIURPRIN2017 project Coarsegrained description for nonequilibrium systems and transport phenomena (CONEST) n.201798CZL. S.J.Y. acknowledges support by the Midcareer Researcher Program (No. 2019R1A2C3006242) through the National Research Foundation of Korea. We thank, A. A. Trani, L. Ciotti, and G. Ciraolo for the discussions at an early stage of this project, and the anonymous Referee for his/her comments that helped improving the presentation of this work.
References
 Aarseth, S. J. 2003, Gravitational NBody Simulations (Cambridge, UK: Cambridge University Press) [Google Scholar]
 Agarwal, M., & Milosavljević, M. 2011, ApJ, 729, 35 [NASA ADS] [CrossRef] [Google Scholar]
 AlfaroCuello, M., Kacharov, N., Neumayer, N., et al. 2020a, in Star Clusters: From the Milky Way to the Early Universe, eds. A. Bragaglia, M. Davies, A. Sills, & E. Vesperini, 351, 47 [NASA ADS] [Google Scholar]
 AlfaroCuello, M., Kacharov, N., Neumayer, N., et al. 2020b, ApJ, 892, 20 [CrossRef] [Google Scholar]
 AmaroSeoane, P., Konstantinidis, S., Brem, P., & Catelan, M. 2013, MNRAS, 435, 809 [NASA ADS] [CrossRef] [Google Scholar]
 Antonini, F., & Gieles, M. 2020, MNRAS, 492, 2936 [CrossRef] [Google Scholar]
 Antonini, F., CapuzzoDolcetta, R., MastrobuonoBattisti, A., & Merritt, D. 2012, ApJ, 750, 111 [NASA ADS] [CrossRef] [Google Scholar]
 ArcaSedda, M., & CapuzzoDolcetta, R. 2014, MNRAS, 444, 3738 [NASA ADS] [CrossRef] [Google Scholar]
 Arca Sedda, M., Askar, A., & Giersz, M. 2018, MNRAS, 479, 4652 [NASA ADS] [CrossRef] [Google Scholar]
 Askar, A., Szkudlarek, M., GondekRosińska, D., Giersz, M., & Bulik, T. 2017, MNRAS, 464, L36 [NASA ADS] [CrossRef] [Google Scholar]
 Askar, A., Davies, M. B., & Church, R. P. 2021, MNRAS, 502, 2682 [NASA ADS] [CrossRef] [Google Scholar]
 Bae, Y.B., Kim, C., & Lee, H. M. 2014, MNRAS, 440, 2714 [NASA ADS] [CrossRef] [Google Scholar]
 Bailyn, C. D. 1995, ARA&A, 33, 133 [NASA ADS] [CrossRef] [Google Scholar]
 Banerjee, S., Baumgardt, H., & Kroupa, P. 2010, MNRAS, 402, 371 [NASA ADS] [CrossRef] [Google Scholar]
 Belczynski, K., Kalogera, V., & Bulik, T. 2002, ApJ, 572, 407 [NASA ADS] [CrossRef] [Google Scholar]
 Belokurov, V., Erkal, D., Evans, N. W., Koposov, S. E., & Deason, A. J. 2018, MNRAS, 478, 611 [Google Scholar]
 Bouffanais, Y., Mapelli, M., Gerosa, D., et al. 2019, ApJ, 886, 25 [NASA ADS] [CrossRef] [Google Scholar]
 Breen, P. G., Foley, C. N., Boekholt, T., & Portegies Zwart, S. 2020, MNRAS, 494, 2465 [Google Scholar]
 Breivik, K., Rodriguez, C. L., Larson, S. L., Kalogera, V., & Rasio, F. A. 2016, ApJ, 830, L18 [NASA ADS] [CrossRef] [Google Scholar]
 CapuzzoDolcetta, R. 1993, ApJ, 415, 616 [NASA ADS] [CrossRef] [Google Scholar]
 CapuzzoDolcetta, R., & Miocchi, P. 2008, MNRAS, 388, L69 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Chatterjee, P., Hernquist, L., & Loeb, A. 2002a, Phys. Rev. Lett., 88, 121103 [CrossRef] [Google Scholar]
 Chatterjee, P., Hernquist, L., & Loeb, A. 2002b, ApJ, 572, 371 [NASA ADS] [CrossRef] [Google Scholar]
 Chatterjee, P., Hernquist, L., & Loeb, A. 2003, ApJ, 592, 32 [CrossRef] [Google Scholar]
 Chatterjee, S., Rodriguez, C. L., Kalogera, V., & Rasio, F. A. 2017, ApJ, 836, L26 [NASA ADS] [CrossRef] [Google Scholar]
 Chung, C., Pasquato, M., Lee, S.Y., et al. 2019, ApJ, 883, L31 [NASA ADS] [CrossRef] [Google Scholar]
 Ciotti, L. 2010, in American Institute of Physics Conference Series, eds. G. Bertin, F. de Luca, G. Lodato, R. Pozzoli, & M. Romé, 1242, 117 [NASA ADS] [Google Scholar]
 Ciotti, L., Nipoti, C., & Londrillo, P. 2007, Proceedings of the International Workshop Collective Phenomena in Macroscopic Systems (Singapore: World Scientific), 177 [Google Scholar]
 Ciraolo, G., Bufferand, H., Di Cintio, P., et al. 2018, Contrib. Plasma Phys., 58, 457 [NASA ADS] [CrossRef] [Google Scholar]
 Cohn, H. 1980, ApJ, 242, 765 [NASA ADS] [CrossRef] [Google Scholar]
 Cole, D. R., Debattista, V. P., Varri, A.L., Adam, M., & Seth, A. C. 2017, MNRAS, 466, 2895 [NASA ADS] [CrossRef] [Google Scholar]
 Di Carlo, U. N., Giacobbo, N., Mapelli, M., et al. 2019, MNRAS, 487, 2947 [NASA ADS] [CrossRef] [Google Scholar]
 Di Cintio, P., & Casetti, L. 2019, MNRAS, 489, 5876 [NASA ADS] [CrossRef] [Google Scholar]
 Di Cintio, P., & Casetti, L. 2020, MNRAS, 494, 1027 [NASA ADS] [CrossRef] [Google Scholar]
 Di Cintio, P., Saalmann, U., & Rost, J.M. 2013, Phys. Rev. Lett., 111, 123401 [NASA ADS] [CrossRef] [Google Scholar]
 Di Cintio, P., Livi, R., Bufferand, H., et al. 2015, Phys. Rev. E, 92, 062108 [NASA ADS] [CrossRef] [Google Scholar]
 Di Cintio, P., Livi, R., Lepri, S., & Ciraolo, G. 2017a, Phys. Rev. E, 95, 043203 [NASA ADS] [CrossRef] [Google Scholar]
 Di Cintio, P., Gupta, S., & Casetti, L. 2017b, Mem. Soc. Astron. It., 88, 733 [NASA ADS] [Google Scholar]
 Di Cintio, P., Gupta, S., & Casetti, L. 2018, MNRAS, 475, 1137 [NASA ADS] [CrossRef] [Google Scholar]
 Di Cintio, P., Ciotti, L., & Nipoti, C. 2020, in Star Clusters: From the Milky Way to the Early Universe, eds. A. Bragaglia, M. Davies, A. Sills, & E. Vesperini, 351, 93 [NASA ADS] [Google Scholar]
 Di Matteo, P., Haywood, M., Lehnert, M. D., et al. 2019, A&A, 632, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ebisuzaki, T., Makino, J., Tsuru, T. G., et al. 2001, ApJ, 562, L19 [NASA ADS] [CrossRef] [Google Scholar]
 Fabian, A. C., Pringle, J. E., & Rees, M. J. 1975, MNRAS, 172, 15 [NASA ADS] [CrossRef] [Google Scholar]
 Fregeau, J. 2012, Astrophysics Source Code Library [record ascl:1208.011] [Google Scholar]
 Fregeau, J. M., Cheung, P., Portegies Zwart, S. F., & Rasio, F. A. 2004, MNRAS, 352, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Freitag, M., & Benz, W. 2001, A&A, 375, 711 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Freitag, M., & Benz, W. 2002, A&A, 394, 345 [NASA ADS] [CrossRef] [EDP Sciences] [PubMed] [Google Scholar]
 Fukushige, T., & Heggie, D. C. 2000, MNRAS, 318, 753 [NASA ADS] [CrossRef] [Google Scholar]
 Giacobbo, N., Mapelli, M., & Spera, M. 2018, MNRAS, 474, 2959 [NASA ADS] [CrossRef] [Google Scholar]
 Giersz, M. 2001, MNRAS, 324, 218 [NASA ADS] [CrossRef] [Google Scholar]
 Giersz, M. 2006, MNRAS, 371, 484 [NASA ADS] [CrossRef] [Google Scholar]
 Giersz, M., Heggie, D. C., Hurley, J. R., & Hypki, A. 2013, MNRAS, 431, 2184 [NASA ADS] [CrossRef] [Google Scholar]
 Gompper, G., Ihle, T., Kroll, D. M., & Winkler, R. G. 2009, in MultiParticle Collision Dynamics: A ParticleBased Mesoscale Simulation Approach to the Hydrodynamics of Complex Fluids, eds. C. Holm, & K. Kremer, 1 [Google Scholar]
 Grand, R. J. J., Kawata, D., Belokurov, V., et al. 2020, MNRAS, 497, 1603 [CrossRef] [Google Scholar]
 Heggie, D. C. 2016, Mem. Soc. Astron. It., 87, 579 [NASA ADS] [Google Scholar]
 Heggie, D. C., & Stevenson, D. 1988, MNRAS, 230, 223 [NASA ADS] [CrossRef] [Google Scholar]
 Hénon, M. 1971a, Ap&SS, 13, 284 [NASA ADS] [CrossRef] [Google Scholar]
 Hénon, M. H. 1971b, Ap&SS, 14, 151 [NASA ADS] [CrossRef] [Google Scholar]
 Hénon, M. 1975, in Dynamics of the Solar Systems, ed. A. Hayli, IAU Symp., 69, 133 [Google Scholar]
 Hockney, R. W., & Eastwood, J. W. 1988, Computer Simulation using Particles (Bristol: Hilger) [Google Scholar]
 Hurley, J. R., & Shara, M. M. 2012, MNRAS, 425, 2872 [NASA ADS] [CrossRef] [Google Scholar]
 Hurley, J. R., Pols, O. R., & Tout, C. A. 2000, MNRAS, 315, 543 [NASA ADS] [CrossRef] [Google Scholar]
 Hurley, J. R., Tout, C. A., & Pols, O. R. 2002, MNRAS, 329, 897 [NASA ADS] [CrossRef] [Google Scholar]
 Hurley, J. R., Sippel, A. C., Tout, C. A., & Aarseth, S. J. 2016, PASA, 33, e036 [NASA ADS] [CrossRef] [Google Scholar]
 Hypki, A., & Giersz, M. 2013, MNRAS, 429, 1221 [NASA ADS] [CrossRef] [Google Scholar]
 Ibata, R. A., Bellazzini, M., Malhan, K., Martin, N., & Bianchini, P. 2019, Nat. Astron., 3, 667 [NASA ADS] [CrossRef] [Google Scholar]
 Joshi, K. J., Rasio, F. A., & Portegies Zwart, S. 2000, ApJ, 540, 969 [NASA ADS] [CrossRef] [Google Scholar]
 Kremer, K., Chatterjee, S., Breivik, K., et al. 2018, Phys. Rev. Lett., 120, 191103 [NASA ADS] [CrossRef] [Google Scholar]
 Kroupa, P. 2002, Science, 295, 82 [Google Scholar]
 Küpper, A. H. W., Kroupa, P., & Baumgardt, H. 2008, MNRAS, 389, 889 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, Y. W., Joo, J. M., Sohn, Y. J., et al. 1999, Nature, 402, 55 [Google Scholar]
 Leigh, N., Sills, A., & Knigge, C. 2007, ApJ, 661, 210 [NASA ADS] [CrossRef] [Google Scholar]
 Lepri, S., Bufferand, H., Ciraolo, G., et al. 2019, in Stochastic Dynamics Out of Equilibrium, eds. G. Giacomin, S. Olla, E. Saada, H. Spohn, & G. Stoltz (Cham: Springer International Publishing), 364 [Google Scholar]
 Makino, J., & Hut, P. 1988, ApJS, 68, 833 [NASA ADS] [CrossRef] [Google Scholar]
 Malevanets, A., & Kapral, R. 1999, J. Chem. Phys., 110, 8605 [CrossRef] [Google Scholar]
 Malevanets, A., & Kapral, R. 2004, in Novel Methods in Soft Matter Simulations, eds. M. Karttunen, A. Lukkarinen, & I. Vattulainen (Berlin: Springer Verlag), Lect. Notes Phys., 640, 116 [NASA ADS] [CrossRef] [Google Scholar]
 Mapelli, M. 2016, MNRAS, 459, 3432 [NASA ADS] [CrossRef] [Google Scholar]
 Massari, D., Koppelman, H. H., & Helmi, A. 2019, A&A, 630, L4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 MastrobuonoBattisti, A., Perets, H. B., & Loeb, A. 2014, ApJ, 796, 40 [NASA ADS] [CrossRef] [Google Scholar]
 McLachlan, R. I., & Atela, P. 1992, Nonlinearity, 5, 541 [CrossRef] [Google Scholar]
 McLaughlin, D. E., & van der Marel, R. P. 2005, ApJS, 161, 304 [NASA ADS] [CrossRef] [Google Scholar]
 Mikkola, S. 2008, in Regular Algorithms for the FewBody Problem, eds. S. J. Aarseth, C. A. Tout, & R. A. Mardling, 760, 31 [NASA ADS] [Google Scholar]
 Miocchi, P., Pasquato, M., Lanzoni, B., et al. 2015, ApJ, 799, 44 [CrossRef] [Google Scholar]
 Misgeld, I., & Hilker, M. 2011, MNRAS, 414, 3699 [Google Scholar]
 Myeong, G. C., Evans, N. W., Belokurov, V., Sanders, J. L., & Koposov, S. E. 2018, MNRAS, 478, 5449 [NASA ADS] [CrossRef] [Google Scholar]
 Myeong, G. C., Vasiliev, E., Iorio, G., Evans, N. W., & Belokurov, V. 2019, MNRAS, 488, 1235 [Google Scholar]
 Neumayer, N. 2017, in Formation, Evolution, and Survival of Massive Star Clusters, eds. C. Charbonnel, & A. Nota, IAU Symp., 316, 84 [Google Scholar]
 Nitadori, K., & Aarseth, S. J. 2012, MNRAS, 424, 545 [NASA ADS] [CrossRef] [Google Scholar]
 Pasquato, M., & Chung, C. 2016, A&A, 589, A95 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pasquato, M., & Di Cintio, P. 2020, A&A, 640, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pasquato, M., de Luca, A., Raimondo, G., et al. 2014, ApJ, 789, 28 [NASA ADS] [CrossRef] [Google Scholar]
 Pasquato, M., Miocchi, P., & Yoon, S.J. 2018, ApJ, 867, 163 [CrossRef] [Google Scholar]
 Pattabiraman, B., Umbreit, S., Liao, W.K., et al. 2013, ApJS, 204, 15 [NASA ADS] [CrossRef] [Google Scholar]
 Perets, H. B., & MastrobuonoBattisti, A. 2014, ApJ, 784, L44 [NASA ADS] [CrossRef] [Google Scholar]
 Pijloo, J. T., Portegies Zwart, S. F., Alexander, P. E. R., et al. 2015, MNRAS, 453, 605 [NASA ADS] [CrossRef] [Google Scholar]
 Plummer, H. C. 1911, MNRAS, 71, 460 [CrossRef] [Google Scholar]
 Portegies Zwart, S. F., McMillan, S. L. W., Hut, P., & Makino, J. 2001, MNRAS, 321, 199 [NASA ADS] [CrossRef] [Google Scholar]
 Portegies Zwart, S. F., Baumgardt, H., Hut, P., Makino, J., & McMillan, S. L. W. 2004, Nature, 428, 724 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Portegies Zwart, S. F., Baumgardt, H., McMillan, S. L. W., et al. 2006, ApJ, 641, 319 [NASA ADS] [CrossRef] [Google Scholar]
 Portegies Zwart, S. F., McMillan, S. L. W., & Gieles, M. 2010, ARA&A, 48, 431 [NASA ADS] [CrossRef] [Google Scholar]
 Rastello, S., AmaroSeoane, P., ArcaSedda, M., et al. 2019, MNRAS, 483, 1233 [NASA ADS] [CrossRef] [Google Scholar]
 Rastello, S., Mapelli, M., Di Carlo, U. N., et al. 2020, MNRAS, 497, 1563 [CrossRef] [Google Scholar]
 Rodriguez, C. L., Chatterjee, S., & Rasio, F. A. 2016, Phys. Rev. D, 93, 084029 [NASA ADS] [CrossRef] [Google Scholar]
 Rodriguez, C. L., Pattabiraman, B., Chatterjee, S., et al. 2018, Comput. Astrophys. Cosmol., 5, 5 [NASA ADS] [CrossRef] [Google Scholar]
 Ryder, J. 2005, PhD Thesis, Oxford University, UK [Google Scholar]
 Salpeter, E. E. 1955, ApJ, 121, 161 [Google Scholar]
 Sarajedini, A., & Layden, A. C. 1995, AJ, 109, 1086 [Google Scholar]
 Sollima, A., & Ferraro, F. R. 2019, MNRAS, 483, 1523 [Google Scholar]
 Sollima, A., & Mastrobuono Battisti, A. 2014, MNRAS, 443, 3513 [NASA ADS] [CrossRef] [Google Scholar]
 Spera, M., & Mapelli, M. 2017, MNRAS, 470, 4739 [Google Scholar]
 Stodolkiewicz, J. S. 1982, Acta Astron., 32, 63 [NASA ADS] [Google Scholar]
 Stodolkiewicz, J. S. 1986, Acta Astron., 36, 19 [NASA ADS] [Google Scholar]
 Takahashi, K., & Portegies Zwart, S. F. 2000, ApJ, 535, 759 [NASA ADS] [CrossRef] [Google Scholar]
 van den Berg, M. 2019, Proc. IAU, 14, 367 [Google Scholar]
 Vasiliev, E. 2015, MNRAS, 446, 3150 [NASA ADS] [CrossRef] [Google Scholar]
 Verbunt, F., & Lewin, W. H. G. 2006, Globular Cluster Xray Sources (Cambridge, UK: Cambridge University Press), 341 [Google Scholar]
 Walcher, C. J., van der Marel, R. P., McLaughlin, D., et al. 2005, ApJ, 618, 237 [NASA ADS] [CrossRef] [Google Scholar]
 Wang, L., Spurzem, R., Aarseth, S., et al. 2016, MNRAS, 458, 1450 [NASA ADS] [CrossRef] [Google Scholar]
 Wang, L., Kroupa, P., Takahashi, K., & Jerabkova, T. 2020, MNRAS, 491, 440 [CrossRef] [Google Scholar]
 Ziosi, B. M., Mapelli, M., Branchesi, M., & Tormen, G. 2014, MNRAS, 441, 3703 [NASA ADS] [CrossRef] [Google Scholar]
All Figures
Fig. 1. Globular star clusters of the Milky Way and satellites from McLaughlin & van der Marel (2005). Log halfmass relaxation time (in years) is plotted against log total mass (in Solar mass). Globular clusters Ω Centauri and M 54 are shown as a magenta square and cyan triangle, respectively. Assuming a mean stellar mass of 0.5 M_{⊙}, star clusters with a mass above 10^{5} M_{⊙} (shaded in light blue) contain 2 × 10^{5} stars and can be simulated only with great computational effort, and they are limited to a handful of realisations. This is especially true if a realistic binary fraction is included, at variance with smaller systems (grey shaded area) that are well within the capabilities of current direct Nbody codes. Star clusters with a relaxation time over the typical globular cluster age (≈10 Gyr) can be regarded as collisionless and are shaded in peach and dark green. According to this definition, Ω Centauri and M 54 are only slightly collisional, but clearly not accessible to modelling through direct Nbody simulations. 

In the text 
Fig. 2. Evolution of the the virial ratio −2K/U (top panel), the fluctuations of the total energy δE in units of the initial energy E_{0} (middle panel), and the norm of the total angular momentum δL (bottom panel) for the same isotropic Plummer initial conditions with N = 2 × 10^{4} particles, evolved with the direct Nbody code (black lines), the exact angular momentum preserving MPC scheme (light blue lines), and the faster MPC rotation with a random axis (purple lines). The initial conditions have a slight angular momentum due to the randomised initialisation procedure for stellar velocities. 

In the text 
Fig. 3. Projections on the x − y plane of two tracer orbits starting from the same initial conditions propagated in a Plummer system with N = 2 × 10^{4} equal mass particles, evolved with an Nbody code (panels a and b, black lines) and MPCDSS (panels c and d, red lines). Bottom panels (e–h): corresponding phasespace sections in the r − v_{r} subspace. 

In the text 
Fig. 4. Fourier spectra of the radial coordinate r for typical orbits extracted from an Nbody simulation (red curve) and an MPC simulation (black curve) for a system of N = 2 × 10^{4} equal mass particles. 

In the text 
Fig. 5. Probability density function of the x coordinate 𝒫(x) for a tracer particle starting with the same initial condition in a N = 2 × 10^{4} Plummer model evolved with the direct Nbody (black circles) and the MPC (red squares) methods. 

In the text 
Fig. 6. Number of escapers N_{esc} as a function of time in units of the dynamical time t_{dyn} for a Plummer model with N = 32 000 and powerlaw mass spectrum with, from left to right, α = 1.5, 2.0, 2.3, 2.5, and 3.0 when evolved with NBODY6 (dashed lines) and with MPCDSS (solid lines). 

In the text 
Fig. 7. Evolution of the central density ρ_{0} (upper panel) and central velocity dispersion σ_{0} (lower panel) for the Nbody simulation (empty symbols) and the MPC simulation (filled symbols) for the model with α = 2.3 in Fig. 6. The symbol sizes are of the order of the mean (Poissonian) error of particle counts in both plots. 

In the text 
Fig. 8. Time of core collapse t_{cc} in units of the dynamical time t_{dyn} as a function of the mass spectrum exponent α for MPC (empty red squares) and Nbody (small black circles) simulations. A parabolic fit to the Nbody simulations (thick green line) is superimposed. The MPC result is generally within the range spanned by the direct Nbody realisations, except for a few values of α at the high end, where core collapse happens earlier in direct Nbody simulations. 

In the text 
Fig. 9. Relations of the number of particles N versus the time of core collapse t_{cc} in units of the dynamical time (left panel) and relaxation time (right panel) for singlecomponent Plummer models. 

In the text 
Fig. 10. 3D density profile for a model with N = 10^{5} equal mass particles and initial Plummer density distribution (thin black solid line) at t = 1, 2, 3, 6, 8, 10, and 18 twobody relaxation times t_{2b} (coloured lines). The heavy dashed line marks the theoretical r^{−2.23} profile. 

In the text 
Fig. 11. Evolution of the 3D Lagrangian radii enclosing 2%, 5%, 10%, 50%, and 90% of the total number of particles N for the same model in Fig. 10. The vertical dashed line marks the system’s corecollapse time t_{cc} ≈ 10 t_{2b} ≈ 12 000 t_{dyn}. 

In the text 
Fig. 12. Evolution of the 3D density profile for two models with N = 10^{5} and initial Plummer density distribution (thin solid line) and mass spectra with α = 1.5 (left panel) and 2.0 (right panel). As in Fig. 10, the heavy dashed line marks the theoretical r^{−2.23} profile. 

In the text 
Fig. 13. Evolution of the 3D Lagrangian radii enclosing 2%, 5%, 10%, 50%, and 90% of the total number of particles N = 10^{5} for the two models shown in Fig. 12. The solid and dotteddashed lines refer the α = 1.5 and 2 cases, respectively. The two vertical dashed lines mark the corecollapse times for the models: t_{cc} ≈ 2000 t_{dyn} for α = 1.5 and t_{cc} ≈ 773 t_{dyn} for α = 2.0. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.