EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 649 (May 2021) A&A, 649 (2021) A24 Full HTML
Free Access
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/0004-6361/202038784
Published online 05 May 2021

© 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), multi-messenger 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 intermediate-mass black hole seeds and a delivery mechanism to the galactic centres; see e.g., Capuzzo-Dolcetta 1993; Ebisuzaki et al. 2001; Portegies Zwart et al. 2004, 2006; Capuzzo-Dolcetta & Miocchi 2008; Antonini et al. 2012; Mastrobuono-Battisti et al. 2014; Arca-Sedda & Capuzzo-Dolcetta 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 self-gravitating N-body systems with a realistic number of stellar particles (in some cases well above 106) over several relaxation times is extremely challenging in terms of computational resources. This is due to the super-quadratic scaling of complexity with the number of particles in direct summation codes (Makino & Hut 1988; Aarseth 2003). Current state-of-the-art direct N-body simulations with 106 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 two-body 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 N-body simulations infeasible, especially when the need to obtain a significant number of realisations of the same system is taken into account.

thumbnail Fig. 1.

Globular star clusters of the Milky Way and satellites from McLaughlin & van der Marel (2005). Log half-mass 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 105 M⊙ (shaded in light blue) contain 2 × 105 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 N-body 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 N-body simulations.

Several approximated alternatives to the direct N-body approach1 that do not share its prohibitive computational cost exist (see Heggie 2016 for an excellent review). The family of so-called Monte Carlo codes, which essentially solve the Fokker-Planck 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 N-body and Monte Carlo (or any Fokker-Planck 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 N-body simulations with N ≈ 2 × 104 when applied to GCs orbiting in a fixed point-like galactic potential.

In this first paper of a series we introduce a new simulation scheme, the Multi-Particle 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 Way2: M 54 and Ω Centauri. Both clusters show a spread in metallicity (Sarajedini & Layden 1995; Lee et al. 1999), hinting at a non-trivial dynamical history that possibly includes one or more mergers (see e.g., Cole et al. 2017; Alfaro-Cuello et al. 2020a,b) that may still affect present-day observable properties (Amaro-Seoane et al. 2013; Pasquato & Chung 2016). The clusters have large masses, thus going beyond the limit of what is currently modellable with ‘honest’ direct N-body 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 state-of-the-art direct N-body 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 multi-particle collision method

In our numerical code, we resolve the collisional interactions between stars using the so-called multi-particle 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 Galilean-invariant dynamics, that the Navier-Stokes 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 non-collisional 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 cell3 onto which the simulation domain has been partitioned.

At the beginning of the collision step, the code evaluates the centre-of-mass (c.o.m.) velocity in every cell:

u com , i = 1 m tot , i ∑ j = 1 n i m j v j ; m tot , i = ∑ j = 1 n i m j , $$ \begin{aligned} {\boldsymbol{u}}_{\mathrm{com},i}=\frac{1}{m_{\mathrm{tot},i}}\sum _{j=1}^{n_i}m_j{\boldsymbol{v}}_j;\quad m_{\mathrm{tot},i}=\sum _{j=1}^{n_i}m_j, \end{aligned} $$(1)

and the relative velocities δvj = vj − ui. Then, for each cell a random axis Ri and rotation angle αi are sampled from uniform distributions. At this point, the vectors δvj are rotated around Ri of αi and then converted back to the simulation frame, so that for the jth particle in cell i the new velocity reads

v j ′ = u i + δ v j , ⊥ cos ( α i ) + ( δ v j , ⊥ × R i ) sin ( α i ) + δ v j , ‖ , $$ \begin{aligned} {\boldsymbol{v}}_{j}^{\prime }={\boldsymbol{u}}_i+\delta {\boldsymbol{v}}_{j,\perp }\mathrm{cos}(\alpha _i)+(\delta {\boldsymbol{v}}_{j,\perp }\times {\boldsymbol{R}}_i)\mathrm{sin}(\alpha _i)+\delta {\boldsymbol{v}}_{j,\parallel }, \end{aligned} $$(2)

where δvj, ⊥ and δvj, ∥ are the relative velocity components perpendicular and parallel to Ri, respectively.

Such an operation conserves exactly the total kinetic energy Ki and the three components of the momentum Pi 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 Li by defining αi such that

sin ( α i ) = − 2 a i b i a i 2 + b i 2 ; cos ( α i ) = a i 2 − b i 2 a i 2 + b i 2 , $$ \begin{aligned} \mathrm{sin}(\alpha _i)=-\frac{2a_ib_i}{a_i^2+b_i^2};\quad \mathrm{cos}(\alpha _i)=\frac{a_i^2-b_i^2}{a_i^2+b_i^2}, \end{aligned} $$(3)

where

a i = ∑ j = 1 N i [ r j × ( v j − u i ) ] | z ; b i = ∑ j = 1 N i r j · ( v j − u i ) . $$ \begin{aligned} a_i=\sum _{j=1}^{N_i}\left[{\boldsymbol{r}}_j\times ({\boldsymbol{v}}_j-{\boldsymbol{u}}_i)\right]|_z;\quad b_i=\sum _{j=1}^{N_i}{\boldsymbol{r}}_j\cdot ({\boldsymbol{v}}_j-{\boldsymbol{u}}_i). \end{aligned} $$(4)

In the formulae above, rj 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

v j ′ = u i + G α i , i · δ v j , $$ \begin{aligned} {\boldsymbol{v}}_{j}^{\prime }={\boldsymbol{u}}_i+{\mathcal{G} }_{\alpha _i,i}\cdot \delta {\boldsymbol{v}}_{j}, \end{aligned} $$(5)

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 Ki and Pi. We also note that the conservation of the total angular momentum can be achieved even in 3D systems by choosing Ri to be parallel to the direction of the cell’s angular momentum vector Li and taking a i = ∑ j = 1 N i [ r j × ( v j − u i ) ] | L i $ a_i=\sum\nolimits_{j=1}^{N_i}\left[\boldsymbol{r}_j\times(\boldsymbol{v}_j-\boldsymbol{u}_i)\right]|_{\boldsymbol{L}_i} $ 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 time-consuming 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 cell-dependent probability accounting for the local degree of collisionality. We first defined the cell-dependent MPC probability as

p i = Erf ( β Δ t 8 π G 2 m ¯ i 2 n ¯ log Λ i σ i 3 ) , $$ \begin{aligned} p_i=\mathrm{Erf}\left(\beta \frac{\Delta t8\pi G^2\bar{m}^2_i\bar{n}\log \Lambda _i}{\sigma ^3_i}\right), \end{aligned} $$(6)

where Δt is the timestep, n ¯ $ \bar{n} $ the mean stellar number density, m ¯ i $ \bar{m}_i $ and σi the average particle mass and the velocity dispersion in the cell, respectively, β is a dimensionless constant fixed to 2Nc, and Erf(x) the standard error function. The cell-dependent Coulomb logarithm is defined as

log Λ i = log ( σ i 2 r s / 2 G m ¯ i ) , $$ \begin{aligned} \log \Lambda _i=\log (\sigma ^2_ir_{\rm s}/2G\bar{m}_i), \end{aligned} $$(7)

with rs 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 multi-particle collision is applied for all cells for which p*i ≤ pi. 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 cell-crossing 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. Vice-versa, 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 particle-mesh scheme (see e.g., Hockney & Eastwood 1988) that solves the Poisson equation:

Δ Φ r = − 4 π G ρ r $$ \begin{aligned} \Delta \Phi _{\boldsymbol{r}}=-4\pi G\rho _{\boldsymbol{r}} \end{aligned} $$(8)

on a spherical grid in polar coordinates Nr × Nϑ × Nφ and interpolates ∇Φ at each particle position ri.

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 tdyn (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

r ¨ i = − G M ( r i ) r i 3 r i , $$ \begin{aligned} \ddot{\boldsymbol{r}}_i=-\frac{GM(r_i)}{r_i^3}{\boldsymbol{r}}_i, \end{aligned} $$(9)

where M(ri) is the mass within the particle’s radial coordinate ri. By doing so, when needed, the potential Φ(ri) on particle i can be obtained as

Φ ( r i ) = − G ( M ( r i ) r i + ∑ j = i + 1 N m j r j ) , $$ \begin{aligned} \Phi (r_i)=-G\left(\frac{M(r_i)}{r_i}+\sum _{j=i+1}^N \frac{m_j}{r_j}\right), \end{aligned} $$(10)

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 Ji = miri × vi would be conserved individually.

3. Comparison with direct N-body

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 N-body 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 N-body 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 − E0) and δL = (L − L0) (where E0 and L0 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:

ρ ( r ) = 3 4 π M r s 2 ( r s 2 + r 2 ) 5 / 2 , $$ \begin{aligned} \rho (r)=\frac{3}{4\pi }\frac{Mr_{\rm s}^2}{(r_{\rm s}^2+r^2)^{5/2}}, \end{aligned} $$(11)

with total mass M and scale radius rs, related to the half-mass radius by rhalf ≈ 1.3 rs. In both the N-body and MPC cases, the systems were evolved up to 1000 dynamical times, defined by

t dyn = r s 3 GM , $$ \begin{aligned} t_{\rm dyn}=\sqrt{\frac{r_{\rm s}^3}{GM}}, \end{aligned} $$(12)

with a fixed timestep Δt = tdyn/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 × 104 particles evolved with the direct N-body code and the two versions of MPCDSS with and without the angular momentum conserving scheme. We only carried out the simulation up to 500 tdyn so that the collisions would not affect the model’s properties much. In all cases, we used the same fixed timestep Δt and a second-order propagation scheme.

thumbnail 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 E0 (middle panel), and the norm of the total angular momentum δL (bottom panel) for the same isotropic Plummer initial conditions with N = 2 × 104 particles, evolved with the direct N-body 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 N-body integration. This is due to the fact that the enforced spherical symmetry and the smoother grid-based 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/E0 as small as 10−3 for the MPC simulations at variance with the direct N-body 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 set-up, its fluctuations δL/L0 are of the order of 10−4, while they are of the order of 10−3 for the direct N-body 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/L0 of the order of 0.5 the smaller the initial value of the angular momentum L0. 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 N-body 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 × 104 integrated with a direct N-body code (panels a and b, black lines) and with MPCDSS (panels c and d, red lines). Always confined within less than two half-mass radii, the particles experience several ‘close encounters’ in both cases, thus subjecting them to dramatic changes in orbital inclination and precession frequency. Direct N-body and MPC dynamics results in a large degree of phase-space ‘diffusion’ with particles exploring the whole energetically accessible region as shown in the lower panels of Fig. 3.

thumbnail 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 × 104 equal mass particles, evolved with an N-body code (panels a and b, black lines) and MPCDSS (panels c and d, red lines). Bottom panels (e–h): corresponding phase-space sections in the r − vr sub-space.

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 N-body 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 × 104 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 × 104). 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 tdyn. 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 N-body orbits has also been made for massive tracer particle orbits propagated semi-analitically with Fokker-Planck 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 = 105 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 Fokker-Planck equation has the best agreement with the results of the direct N-body 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.

thumbnail Fig. 4.

Fourier spectra of the radial coordinate r for typical orbits extracted from an N-body simulation (red curve) and an MPC simulation (black curve) for a system of N = 2 × 104 equal mass particles.

thumbnail Fig. 5.

Probability density function of the x coordinate 𝒫(x) for a tracer particle starting with the same initial condition in a N = 2 × 104 Plummer model evolved with the direct N-body (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 tdyn takes a few days on a dedicated GPU workstation for the direct N-body simulation without stellar evolution and adaptive time step. When using an Nc = 32 × 16 × 16 grid in polar coordinates to perform the collisions and an Nr = 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 multi-mass 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 multi-slope Kroupa (2002) mass function, in this work the particle masses mi were extracted from a pure power-law mass function of the form

F ( m ) = C m α ; m min ≤ m ≤ m max , $$ \begin{aligned} \mathcal{F} (m) = \frac{C}{m^{\alpha }}; \quad m_{\rm min} \le m\le m_{\rm max}, \end{aligned} $$(13)

where α > 0, and the normalisation constant C depends on the minimum-to-maximum mass ratio ℛ = mmin/mmax so that ∫ m min m max F (m)dm=M $ \int_{m_{\rm min}}^{m_{\rm max}}\mathcal{F}(m){\rm d}m=M $.

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 N-body 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 N-body and MPC cases, the systems were evolved for 104 dynamical times, corresponding to roughly 20 two-body (collisional) relaxation times of a model with the same total mass and number of equal mass particles, given by

t 2 b ≈ 0.138 N log N t dyn . $$ \begin{aligned} t_{\rm 2b}\approx \frac{0.138N}{\log N}t_{\rm dyn}. \end{aligned} $$(14)

Again, as a rule, in all sets of MPCDSS simulations we fixed the timestep Δt to tdyn/100. We neglected stellar evolution and the contribution of binaries4. For the two simulation sets (N-body and MPC), we extracted and compared the number of escapers Nesc (defined as the number of particles at r > 17 rs 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 Nesc for choices of α in Eq. (13). We find that in all cases the MPC evolution (solid line) recovers the quasi-linear trend of Nesc with time. However, some discrepancies between MPC and N-body simulations are observed, in particular at low values of α (i.e. ‘flatter’ mass spectrum). We interpret such a difference as the effect of strong two-body collisions between light and heavy particles, resolved in a direct N-body code but somewhat smeared-out in a multi-particle 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 inter-particle collisions with a single move acting on all particles of the cell, the contribution of few but strong light-heavy particle encounters is obviously underestimated.

thumbnail Fig. 6.

Number of escapers Nesc as a function of time in units of the dynamical time tdyn for a Plummer model with N = 32 000 and power-law 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.

thumbnail Fig. 7.

Evolution of the central density ρ0 (upper panel) and central velocity dispersion σ0 (lower panel) for the N-body 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 tcc as the time at which the minimum value reached by r2% (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 mass-function power-law 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 core-collapse. Interestingly, for large values of α (say ≳2.6), tcc starts increasing again in both MPC and direct N-body simulations. We speculate that such curious behaviour might be due to the less efficient dynamical friction in models dominated by low-mass particles (see Ciotti 2010) resulting in larger sinking timescales for particles sitting in the large-mass part of the spectrum. To guide the eye, we fitted the relation between α and tcc using a second-order polynomial (solid green line). We performed additional simulations (to be published elsewhere) with different values of N, α, and ℛ, and the non-monotonic trend of tcc 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 N-body 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 tcc.

thumbnail Fig. 8.

Time of core collapse tcc in units of the dynamical time tdyn as a function of the mass spectrum exponent α for MPC (empty red squares) and N-body (small black circles) simulations. A parabolic fit to the N-body simulations (thick green line) is superimposed. The MPC result is generally within the range spanned by the direct N-body realisations, except for a few values of α at the high end, where core collapse happens earlier in direct N-body 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 single-component models starting with isotropic Plummer density profiles in a broad range of system sizes spanning from 103 to 106. 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 tdyn), as shown in the left panel of Fig. 9 for the equal masses case, with N ranging from 104 to 106. Remarkably, for low values of N, there seems to be a somewhat non-monotonic trend of tcc with N. We performed additional MPC simulations with N as small as 104 using different realisations of the initial condition with choices of timestep and grid size, and such a trend persists.

thumbnail Fig. 9.

Relations of the number of particles N versus the time of core collapse tcc in units of the dynamical time (left panel) and relaxation time (right panel) for single-component Plummer models.

When expressing tcc in units of the collisional relaxation time t2b (cf. Eq. (14)), the picture is inverted and large N systems reach core collapse earlier in units of their intrinsic t2b (see right panel of the same figure). As N increases to larger and larger values, the systems start to approach the so-called 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, t2b also increases as N/log N, resulting in a decreasing trend of tcc/t2b with N.

In Fig. 10, we show the radial density profile at different times between 1 t2b and 18 t2b for the N = 105 case. It is evident, as at around 2 t2b (corresponding roughly to 2340 tdyn 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 t2b, the re-expansion 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 = 104 and as large as N = 106, 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.

thumbnail Fig. 10.

3D density profile for a model with N = 105 equal mass particles and initial Plummer density distribution (thin black solid line) at t = 1, 2, 3, 6, 8, 10, and 18 two-body relaxation times t2b (coloured lines). The heavy dashed line marks the theoretical r−2.23 profile.

We find a core-collapse time tcc (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 t2b (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 N-body simulations of Küpper et al. (2008), finding values between 10 and 15 t2b for models with initial conditions analogous to the ones used in our simulations (i.e. Plummer profile, N = 105, and no mass spectrum).

thumbnail 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 core-collapse time tcc ≈ 10 t2b ≈ 12 000 tdyn.

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 = 105.

thumbnail Fig. 12.

Evolution of the 3D density profile for two models with N = 105 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 tdyn 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 tcc is well below 103 tdyn, being roughly 2000 tdyn for α = 1.5 and 773 tdyn for α = 2.0 as marked by the vertical lines in the figure.

thumbnail Fig. 13.

Evolution of the 3D Lagrangian radii enclosing 2%, 5%, 10%, 50%, and 90% of the total number of particles N = 105 for the two models shown in Fig. 12. The solid and dotted-dashed lines refer the α = 1.5 and 2 cases, respectively. The two vertical dashed lines mark the core-collapse times for the models: tcc ≈ 2000 tdyn for α = 1.5 and tcc ≈ 773 tdyn for α = 2.0.

Since the time-dependent 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 core-collapse 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 direct-summation N-body 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 cell-based rotations of the stellar velocity vectors that, by construction, conserve mass, momentum, and energy5. The MPC approach avoids the intricacies of two- and multiple-body encounters that require indirect N-body code techniques such as softening or Kustaanheimo-Stiefel regularisation (see e.g., Mikkola 2008) in order to treat the singularity of the 1/r gravitational potential6 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 long-term 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 rotating7, 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 N-body 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 N-body 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 N-body 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 intermediate-mass black holes in the Galactic NSC.

1

We note that direct N-body 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).

2

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 N-body simulations (Agarwal & Milosavljević 2011; Perets & Mastrobuono-Battisti 2014) in the context of the so-called repeated accretion scenario (Antonini et al. 2012; Arca-Sedda & Capuzzo-Dolcetta 2014).

3

In the original implementation of the hybrid plasma PIC-MPC code, the cell structure is the same as the one used by the Maxwell solver routines to compute electromagnetic fields.

4

We note that direct N-body simulations, even if started without binaries, may form a few new binaries dynamically.

5

In our implementation of the method, angular momentum conservation is also ensured, but this feature can be switched off to speed up calculations if necessary.

6

We note that adaptively softening the gravitational interaction with 1 / ϵ 2 + r 2 $ 1/\sqrt{\epsilon^2+r^2} $ or preforming the Kustaanheimo-Stiefel mapping of position r and time t in u = r $ u=\sqrt{r} $ and dt = rdτ every time a two-body close encounter takes place, further increases computational time in direct simulations.

7

We note that, since the MPC operator acts in each cell’s c.o.m frame, the presence of a collective rotation field does not affect the stochastic collision mechanism in any way.

Acknowledgments

This project has received funding from the European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie 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 MIUR-PRIN2017 project Coarse-grained description for non-equilibrium systems and transport phenomena (CO-NEST) n.201798CZL. S.-J.Y. acknowledges support by the Mid-career 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

  1. Aarseth, S. J. 2003, Gravitational N-Body Simulations (Cambridge, UK: Cambridge University Press) [Google Scholar]
  2. Agarwal, M., & Milosavljević, M. 2011, ApJ, 729, 35 [NASA ADS] [CrossRef] [Google Scholar]
  3. Alfaro-Cuello, 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]
  4. Alfaro-Cuello, M., Kacharov, N., Neumayer, N., et al. 2020b, ApJ, 892, 20 [CrossRef] [Google Scholar]
  5. Amaro-Seoane, P., Konstantinidis, S., Brem, P., & Catelan, M. 2013, MNRAS, 435, 809 [NASA ADS] [CrossRef] [Google Scholar]
  6. Antonini, F., & Gieles, M. 2020, MNRAS, 492, 2936 [CrossRef] [Google Scholar]
  7. Antonini, F., Capuzzo-Dolcetta, R., Mastrobuono-Battisti, A., & Merritt, D. 2012, ApJ, 750, 111 [NASA ADS] [CrossRef] [Google Scholar]
  8. Arca-Sedda, M., & Capuzzo-Dolcetta, R. 2014, MNRAS, 444, 3738 [NASA ADS] [CrossRef] [Google Scholar]
  9. Arca Sedda, M., Askar, A., & Giersz, M. 2018, MNRAS, 479, 4652 [NASA ADS] [CrossRef] [Google Scholar]
  10. Askar, A., Szkudlarek, M., Gondek-Rosińska, D., Giersz, M., & Bulik, T. 2017, MNRAS, 464, L36 [NASA ADS] [CrossRef] [Google Scholar]
  11. Askar, A., Davies, M. B., & Church, R. P. 2021, MNRAS, 502, 2682 [NASA ADS] [CrossRef] [Google Scholar]
  12. Bae, Y.-B., Kim, C., & Lee, H. M. 2014, MNRAS, 440, 2714 [NASA ADS] [CrossRef] [Google Scholar]
  13. Bailyn, C. D. 1995, ARA&A, 33, 133 [NASA ADS] [CrossRef] [Google Scholar]
  14. Banerjee, S., Baumgardt, H., & Kroupa, P. 2010, MNRAS, 402, 371 [NASA ADS] [CrossRef] [Google Scholar]
  15. Belczynski, K., Kalogera, V., & Bulik, T. 2002, ApJ, 572, 407 [NASA ADS] [CrossRef] [Google Scholar]
  16. Belokurov, V., Erkal, D., Evans, N. W., Koposov, S. E., & Deason, A. J. 2018, MNRAS, 478, 611 [Google Scholar]
  17. Bouffanais, Y., Mapelli, M., Gerosa, D., et al. 2019, ApJ, 886, 25 [NASA ADS] [CrossRef] [Google Scholar]
  18. Breen, P. G., Foley, C. N., Boekholt, T., & Portegies Zwart, S. 2020, MNRAS, 494, 2465 [Google Scholar]
  19. Breivik, K., Rodriguez, C. L., Larson, S. L., Kalogera, V., & Rasio, F. A. 2016, ApJ, 830, L18 [NASA ADS] [CrossRef] [Google Scholar]
  20. Capuzzo-Dolcetta, R. 1993, ApJ, 415, 616 [NASA ADS] [CrossRef] [Google Scholar]
  21. Capuzzo-Dolcetta, R., & Miocchi, P. 2008, MNRAS, 388, L69 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  22. Chatterjee, P., Hernquist, L., & Loeb, A. 2002a, Phys. Rev. Lett., 88, 121103 [CrossRef] [Google Scholar]
  23. Chatterjee, P., Hernquist, L., & Loeb, A. 2002b, ApJ, 572, 371 [NASA ADS] [CrossRef] [Google Scholar]
  24. Chatterjee, P., Hernquist, L., & Loeb, A. 2003, ApJ, 592, 32 [CrossRef] [Google Scholar]
  25. Chatterjee, S., Rodriguez, C. L., Kalogera, V., & Rasio, F. A. 2017, ApJ, 836, L26 [NASA ADS] [CrossRef] [Google Scholar]
  26. Chung, C., Pasquato, M., Lee, S.-Y., et al. 2019, ApJ, 883, L31 [NASA ADS] [CrossRef] [Google Scholar]
  27. 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]
  28. Ciotti, L., Nipoti, C., & Londrillo, P. 2007, Proceedings of the International Workshop Collective Phenomena in Macroscopic Systems (Singapore: World Scientific), 177 [Google Scholar]
  29. Ciraolo, G., Bufferand, H., Di Cintio, P., et al. 2018, Contrib. Plasma Phys., 58, 457 [NASA ADS] [CrossRef] [Google Scholar]
  30. Cohn, H. 1980, ApJ, 242, 765 [NASA ADS] [CrossRef] [Google Scholar]
  31. Cole, D. R., Debattista, V. P., Varri, A.-L., Adam, M., & Seth, A. C. 2017, MNRAS, 466, 2895 [NASA ADS] [CrossRef] [Google Scholar]
  32. Di Carlo, U. N., Giacobbo, N., Mapelli, M., et al. 2019, MNRAS, 487, 2947 [NASA ADS] [CrossRef] [Google Scholar]
  33. Di Cintio, P., & Casetti, L. 2019, MNRAS, 489, 5876 [NASA ADS] [CrossRef] [Google Scholar]
  34. Di Cintio, P., & Casetti, L. 2020, MNRAS, 494, 1027 [NASA ADS] [CrossRef] [Google Scholar]
  35. Di Cintio, P., Saalmann, U., & Rost, J.-M. 2013, Phys. Rev. Lett., 111, 123401 [NASA ADS] [CrossRef] [Google Scholar]
  36. Di Cintio, P., Livi, R., Bufferand, H., et al. 2015, Phys. Rev. E, 92, 062108 [NASA ADS] [CrossRef] [Google Scholar]
  37. Di Cintio, P., Livi, R., Lepri, S., & Ciraolo, G. 2017a, Phys. Rev. E, 95, 043203 [NASA ADS] [CrossRef] [Google Scholar]
  38. Di Cintio, P., Gupta, S., & Casetti, L. 2017b, Mem. Soc. Astron. It., 88, 733 [NASA ADS] [Google Scholar]
  39. Di Cintio, P., Gupta, S., & Casetti, L. 2018, MNRAS, 475, 1137 [NASA ADS] [CrossRef] [Google Scholar]
  40. 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]
  41. Di Matteo, P., Haywood, M., Lehnert, M. D., et al. 2019, A&A, 632, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  42. Ebisuzaki, T., Makino, J., Tsuru, T. G., et al. 2001, ApJ, 562, L19 [NASA ADS] [CrossRef] [Google Scholar]
  43. Fabian, A. C., Pringle, J. E., & Rees, M. J. 1975, MNRAS, 172, 15 [NASA ADS] [CrossRef] [Google Scholar]
  44. Fregeau, J. 2012, Astrophysics Source Code Library [record ascl:1208.011] [Google Scholar]
  45. Fregeau, J. M., Cheung, P., Portegies Zwart, S. F., & Rasio, F. A. 2004, MNRAS, 352, 1 [NASA ADS] [CrossRef] [Google Scholar]
  46. Freitag, M., & Benz, W. 2001, A&A, 375, 711 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  47. Freitag, M., & Benz, W. 2002, A&A, 394, 345 [NASA ADS] [CrossRef] [EDP Sciences] [PubMed] [Google Scholar]
  48. Fukushige, T., & Heggie, D. C. 2000, MNRAS, 318, 753 [NASA ADS] [CrossRef] [Google Scholar]
  49. Giacobbo, N., Mapelli, M., & Spera, M. 2018, MNRAS, 474, 2959 [NASA ADS] [CrossRef] [Google Scholar]
  50. Giersz, M. 2001, MNRAS, 324, 218 [NASA ADS] [CrossRef] [Google Scholar]
  51. Giersz, M. 2006, MNRAS, 371, 484 [NASA ADS] [CrossRef] [Google Scholar]
  52. Giersz, M., Heggie, D. C., Hurley, J. R., & Hypki, A. 2013, MNRAS, 431, 2184 [NASA ADS] [CrossRef] [Google Scholar]
  53. Gompper, G., Ihle, T., Kroll, D. M., & Winkler, R. G. 2009, in Multi-Particle Collision Dynamics: A Particle-Based Mesoscale Simulation Approach to the Hydrodynamics of Complex Fluids, eds. C. Holm, & K. Kremer, 1 [Google Scholar]
  54. Grand, R. J. J., Kawata, D., Belokurov, V., et al. 2020, MNRAS, 497, 1603 [CrossRef] [Google Scholar]
  55. Heggie, D. C. 2016, Mem. Soc. Astron. It., 87, 579 [NASA ADS] [Google Scholar]
  56. Heggie, D. C., & Stevenson, D. 1988, MNRAS, 230, 223 [NASA ADS] [CrossRef] [Google Scholar]
  57. Hénon, M. 1971a, Ap&SS, 13, 284 [NASA ADS] [CrossRef] [Google Scholar]
  58. Hénon, M. H. 1971b, Ap&SS, 14, 151 [NASA ADS] [CrossRef] [Google Scholar]
  59. Hénon, M. 1975, in Dynamics of the Solar Systems, ed. A. Hayli, IAU Symp., 69, 133 [Google Scholar]
  60. Hockney, R. W., & Eastwood, J. W. 1988, Computer Simulation using Particles (Bristol: Hilger) [Google Scholar]
  61. Hurley, J. R., & Shara, M. M. 2012, MNRAS, 425, 2872 [NASA ADS] [CrossRef] [Google Scholar]
  62. Hurley, J. R., Pols, O. R., & Tout, C. A. 2000, MNRAS, 315, 543 [NASA ADS] [CrossRef] [Google Scholar]
  63. Hurley, J. R., Tout, C. A., & Pols, O. R. 2002, MNRAS, 329, 897 [NASA ADS] [CrossRef] [Google Scholar]
  64. Hurley, J. R., Sippel, A. C., Tout, C. A., & Aarseth, S. J. 2016, PASA, 33, e036 [NASA ADS] [CrossRef] [Google Scholar]
  65. Hypki, A., & Giersz, M. 2013, MNRAS, 429, 1221 [NASA ADS] [CrossRef] [Google Scholar]
  66. Ibata, R. A., Bellazzini, M., Malhan, K., Martin, N., & Bianchini, P. 2019, Nat. Astron., 3, 667 [NASA ADS] [CrossRef] [Google Scholar]
  67. Joshi, K. J., Rasio, F. A., & Portegies Zwart, S. 2000, ApJ, 540, 969 [NASA ADS] [CrossRef] [Google Scholar]
  68. Kremer, K., Chatterjee, S., Breivik, K., et al. 2018, Phys. Rev. Lett., 120, 191103 [NASA ADS] [CrossRef] [Google Scholar]
  69. Kroupa, P. 2002, Science, 295, 82 [Google Scholar]
  70. Küpper, A. H. W., Kroupa, P., & Baumgardt, H. 2008, MNRAS, 389, 889 [NASA ADS] [CrossRef] [Google Scholar]
  71. Lee, Y. W., Joo, J. M., Sohn, Y. J., et al. 1999, Nature, 402, 55 [Google Scholar]
  72. Leigh, N., Sills, A., & Knigge, C. 2007, ApJ, 661, 210 [NASA ADS] [CrossRef] [Google Scholar]
  73. 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]
  74. Makino, J., & Hut, P. 1988, ApJS, 68, 833 [NASA ADS] [CrossRef] [Google Scholar]
  75. Malevanets, A., & Kapral, R. 1999, J. Chem. Phys., 110, 8605 [CrossRef] [Google Scholar]
  76. 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]
  77. Mapelli, M. 2016, MNRAS, 459, 3432 [NASA ADS] [CrossRef] [Google Scholar]
  78. Massari, D., Koppelman, H. H., & Helmi, A. 2019, A&A, 630, L4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  79. Mastrobuono-Battisti, A., Perets, H. B., & Loeb, A. 2014, ApJ, 796, 40 [NASA ADS] [CrossRef] [Google Scholar]
  80. McLachlan, R. I., & Atela, P. 1992, Nonlinearity, 5, 541 [CrossRef] [Google Scholar]
  81. McLaughlin, D. E., & van der Marel, R. P. 2005, ApJS, 161, 304 [NASA ADS] [CrossRef] [Google Scholar]
  82. Mikkola, S. 2008, in Regular Algorithms for the Few-Body Problem, eds. S. J. Aarseth, C. A. Tout, & R. A. Mardling, 760, 31 [NASA ADS] [Google Scholar]
  83. Miocchi, P., Pasquato, M., Lanzoni, B., et al. 2015, ApJ, 799, 44 [CrossRef] [Google Scholar]
  84. Misgeld, I., & Hilker, M. 2011, MNRAS, 414, 3699 [Google Scholar]
  85. Myeong, G. C., Evans, N. W., Belokurov, V., Sanders, J. L., & Koposov, S. E. 2018, MNRAS, 478, 5449 [NASA ADS] [CrossRef] [Google Scholar]
  86. Myeong, G. C., Vasiliev, E., Iorio, G., Evans, N. W., & Belokurov, V. 2019, MNRAS, 488, 1235 [Google Scholar]
  87. Neumayer, N. 2017, in Formation, Evolution, and Survival of Massive Star Clusters, eds. C. Charbonnel, & A. Nota, IAU Symp., 316, 84 [Google Scholar]
  88. Nitadori, K., & Aarseth, S. J. 2012, MNRAS, 424, 545 [NASA ADS] [CrossRef] [Google Scholar]
  89. Pasquato, M., & Chung, C. 2016, A&A, 589, A95 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  90. Pasquato, M., & Di Cintio, P. 2020, A&A, 640, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  91. Pasquato, M., de Luca, A., Raimondo, G., et al. 2014, ApJ, 789, 28 [NASA ADS] [CrossRef] [Google Scholar]
  92. Pasquato, M., Miocchi, P., & Yoon, S.-J. 2018, ApJ, 867, 163 [CrossRef] [Google Scholar]
  93. Pattabiraman, B., Umbreit, S., Liao, W.-K., et al. 2013, ApJS, 204, 15 [NASA ADS] [CrossRef] [Google Scholar]
  94. Perets, H. B., & Mastrobuono-Battisti, A. 2014, ApJ, 784, L44 [NASA ADS] [CrossRef] [Google Scholar]
  95. Pijloo, J. T., Portegies Zwart, S. F., Alexander, P. E. R., et al. 2015, MNRAS, 453, 605 [NASA ADS] [CrossRef] [Google Scholar]
  96. Plummer, H. C. 1911, MNRAS, 71, 460 [CrossRef] [Google Scholar]
  97. Portegies Zwart, S. F., McMillan, S. L. W., Hut, P., & Makino, J. 2001, MNRAS, 321, 199 [NASA ADS] [CrossRef] [Google Scholar]
  98. Portegies Zwart, S. F., Baumgardt, H., Hut, P., Makino, J., & McMillan, S. L. W. 2004, Nature, 428, 724 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  99. Portegies Zwart, S. F., Baumgardt, H., McMillan, S. L. W., et al. 2006, ApJ, 641, 319 [NASA ADS] [CrossRef] [Google Scholar]
  100. Portegies Zwart, S. F., McMillan, S. L. W., & Gieles, M. 2010, ARA&A, 48, 431 [NASA ADS] [CrossRef] [Google Scholar]
  101. Rastello, S., Amaro-Seoane, P., Arca-Sedda, M., et al. 2019, MNRAS, 483, 1233 [NASA ADS] [CrossRef] [Google Scholar]
  102. Rastello, S., Mapelli, M., Di Carlo, U. N., et al. 2020, MNRAS, 497, 1563 [CrossRef] [Google Scholar]
  103. Rodriguez, C. L., Chatterjee, S., & Rasio, F. A. 2016, Phys. Rev. D, 93, 084029 [NASA ADS] [CrossRef] [Google Scholar]
  104. Rodriguez, C. L., Pattabiraman, B., Chatterjee, S., et al. 2018, Comput. Astrophys. Cosmol., 5, 5 [NASA ADS] [CrossRef] [Google Scholar]
  105. Ryder, J. 2005, PhD Thesis, Oxford University, UK [Google Scholar]
  106. Salpeter, E. E. 1955, ApJ, 121, 161 [Google Scholar]
  107. Sarajedini, A., & Layden, A. C. 1995, AJ, 109, 1086 [Google Scholar]
  108. Sollima, A., & Ferraro, F. R. 2019, MNRAS, 483, 1523 [Google Scholar]
  109. Sollima, A., & Mastrobuono Battisti, A. 2014, MNRAS, 443, 3513 [NASA ADS] [CrossRef] [Google Scholar]
  110. Spera, M., & Mapelli, M. 2017, MNRAS, 470, 4739 [Google Scholar]
  111. Stodolkiewicz, J. S. 1982, Acta Astron., 32, 63 [NASA ADS] [Google Scholar]
  112. Stodolkiewicz, J. S. 1986, Acta Astron., 36, 19 [NASA ADS] [Google Scholar]
  113. Takahashi, K., & Portegies Zwart, S. F. 2000, ApJ, 535, 759 [NASA ADS] [CrossRef] [Google Scholar]
  114. van den Berg, M. 2019, Proc. IAU, 14, 367 [Google Scholar]
  115. Vasiliev, E. 2015, MNRAS, 446, 3150 [NASA ADS] [CrossRef] [Google Scholar]
  116. Verbunt, F., & Lewin, W. H. G. 2006, Globular Cluster X-ray Sources (Cambridge, UK: Cambridge University Press), 341 [Google Scholar]
  117. Walcher, C. J., van der Marel, R. P., McLaughlin, D., et al. 2005, ApJ, 618, 237 [NASA ADS] [CrossRef] [Google Scholar]
  118. Wang, L., Spurzem, R., Aarseth, S., et al. 2016, MNRAS, 458, 1450 [NASA ADS] [CrossRef] [Google Scholar]
  119. Wang, L., Kroupa, P., Takahashi, K., & Jerabkova, T. 2020, MNRAS, 491, 440 [CrossRef] [Google Scholar]
  120. Ziosi, B. M., Mapelli, M., Branchesi, M., & Tormen, G. 2014, MNRAS, 441, 3703 [NASA ADS] [CrossRef] [Google Scholar]

All Figures

thumbnail Fig. 1.

Globular star clusters of the Milky Way and satellites from McLaughlin & van der Marel (2005). Log half-mass 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 105 M⊙ (shaded in light blue) contain 2 × 105 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 N-body 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 N-body simulations.
In the text
thumbnail 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 E0 (middle panel), and the norm of the total angular momentum δL (bottom panel) for the same isotropic Plummer initial conditions with N = 2 × 104 particles, evolved with the direct N-body 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
thumbnail 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 × 104 equal mass particles, evolved with an N-body code (panels a and b, black lines) and MPCDSS (panels c and d, red lines). Bottom panels (e–h): corresponding phase-space sections in the r − vr sub-space.
In the text
thumbnail Fig. 4.

Fourier spectra of the radial coordinate r for typical orbits extracted from an N-body simulation (red curve) and an MPC simulation (black curve) for a system of N = 2 × 104 equal mass particles.
In the text
thumbnail Fig. 5.

Probability density function of the x coordinate 𝒫(x) for a tracer particle starting with the same initial condition in a N = 2 × 104 Plummer model evolved with the direct N-body (black circles) and the MPC (red squares) methods.
In the text
thumbnail Fig. 6.

Number of escapers Nesc as a function of time in units of the dynamical time tdyn for a Plummer model with N = 32 000 and power-law 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
thumbnail Fig. 7.

Evolution of the central density ρ0 (upper panel) and central velocity dispersion σ0 (lower panel) for the N-body 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
thumbnail Fig. 8.

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

Relations of the number of particles N versus the time of core collapse tcc in units of the dynamical time (left panel) and relaxation time (right panel) for single-component Plummer models.
In the text
thumbnail Fig. 10.

3D density profile for a model with N = 105 equal mass particles and initial Plummer density distribution (thin black solid line) at t = 1, 2, 3, 6, 8, 10, and 18 two-body relaxation times t2b (coloured lines). The heavy dashed line marks the theoretical r−2.23 profile.
In the text
thumbnail 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 core-collapse time tcc ≈ 10 t2b ≈ 12 000 tdyn.
In the text
thumbnail Fig. 12.

Evolution of the 3D density profile for two models with N = 105 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
thumbnail Fig. 13.

Evolution of the 3D Lagrangian radii enclosing 2%, 5%, 10%, 50%, and 90% of the total number of particles N = 105 for the two models shown in Fig. 12. The solid and dotted-dashed lines refer the α = 1.5 and 2 cases, respectively. The two vertical dashed lines mark the core-collapse times for the models: tcc ≈ 2000 tdyn for α = 1.5 and tcc ≈ 773 tdyn for α = 2.0.
In the text

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.