Issue 
A&A
Volume 673, May 2023



Article Number  A8  
Number of page(s)  11  
Section  Galactic structure, stellar clusters and populations  
DOI  https://doi.org/10.1051/00046361/202346124  
Published online  24 April 2023 
Dynamics of intermediate mass black holes in globular clusters
Wander radius and anisotropy profiles
^{1}
Institute of Complex Systems – National Council of Research (ISCCNR), via della Lastruccia 10, 50019 Sesto Fiorentino, Italy
email: pierfrancesco.dicintio@cnr.it
^{2}
Physics and Astronomy Department, University of Firenze, via G. Sansone 1, 50019 Sesto Fiorentino, Italy
^{3}
INAF, Osservatorio Astrofisico di Arcetri, Largo Enrico Fermi, 5, 50125 Firenze, Italy
^{4}
INFN – Sezione di Firenze, via G. Sansone 1, 50019 Sesto Fiorentino, Italy
^{5}
Physics and Astronomy Department Galileo Galilei, University of Padova, Vicolo dell’Osservatorio 3, 35122 Padova, Italy
^{6}
Département de Physique, Université de Montréal, Montreal, Quebec H3T 1J4, Canada
^{7}
The University of Tokyo, Earth Science and Astronomy Department, 731 Hongo, Bunkyoku, Tokyo 1130033, Japan
^{8}
Okinawa Institute of Science and Technology Graduate University, 19191 Tancha, Onnason, Kunigamigun, 9040495 Okinawa, Japan
^{9}
McWilliams Center for Cosmology and Department of Physics, Carnegie Mellon University, Pittsburgh, PA 15213, USA
Received:
10
February
2023
Accepted:
24
March
2023
Context. We recently introduced a new method for simulating collisional gravitational Nbody systems with approximately linear time scaling with N. Our method is based on the multiparticle collision (MPC) scheme, previously applied in fluid dynamics and plasma physics. We were able to simulate globular clusters with a realistic number of stellar particles (at least up to several times 10^{6}) on a standard workstation.
Aims. We simulated clusters hosting an intermediate mass black hole (IMBH), probing a broad range of BHcluster and BH–averagestar mass ratios, unrestricted by the computational constraints that affect direct Nbody codes.
Methods. We set up a grid of hybrid particleincellMPC Nbody simulations using our implementation of the MPC method, MPCDSS. We used either single mass models or models with a Salpeter mass function (a single power law with an exponent of −2.35), with the IMBH initially sitting at the centre. The force exerted by and on the IMBH was evaluated with a direct sum scheme with or without softening. For all simulations we measured the evolution of the Lagrangian radii and core density and velocity dispersion over time. In addition, we also measured the evolution of the velocity anisotropy profiles.
Results. We find that models with an IMBH undergo core collapse at earlier times, the larger the IMBH mass the shallower they are, with an approximately constant central density at core collapse. The presence of an IMBH tends to lower the central velocity dispersion. These results hold independently of the mass function of the model. For the models with Salpeter MF, we observed that equipartition of kinetic energies is never achieved, even long after core collapse. Orbital anisotropy at large radii appears to be driven by energetic escapers on radial orbits, triggered by strong collisions with the IMBH in the core. We measured the wander radius, that is the distance of the IMBH from the centre of mass of the parent system over time, finding that its distribution has positive kurtosis.
Conclusions. Among the results we obtained, which mostly confirm or extend previously known trends that had been established over the range of parameters accessible to direct Nbody simulations, we underline that the leptokurtic nature of the IMBH wander radius distribution might lead to IMBHs presenting as offcentre more frequently than expected, with implications on observational IMBH detection.
Key words: galaxies: star clusters: general / methods: numerical / black hole physics / stars: black holes
© The Authors 2023
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
1. Introduction
A plausible mechanism for super massive black hole (SMBH) seeding is required to explain the observation of quasars at high redshift (Inayoshi et al. 2020; Pacucci & Loeb 2022). Early seeding would rely on pristine gas and is speculated to take place either through direct collapse (Loeb & Rasio 1994; Lodato & Natarajan 2006; Volonteri et al. 2008) or population III stars (Carr et al. 1984; Yoshida et al. 2006; Greif 2015). Later or continuous seeding would instead happen through dynamically mediated gravitational runaway scenarios in dense environments (Miller & Hamilton 2002; Ebisuzaki et al. 2001; Portegies Zwart et al. 2004; Li 2022). SMBH seeds should be detectable today as intermediate mass black holes (IMBHs) in dense stellar systems such as star clusters (Greene et al. 2020; Di Carlo et al. 2021; Rizzuto et al. 2021, 2023), especially if the second scenario is prevalent, modulo expulsion from the host system via gravitational wave recoil kicks (see e.g., HolleyBockelmann et al. 2008; Weatherford et al. 2023). Quantitatively addressing the seeding mechanism requires us to constrain the fraction of ‘rogue’ IMBHs (i.e. not associated with a host star cluster, which may still be detectable by other means, e.g., Ballone et al. 2018), requiring IMBH ejection from the parent cluster to be modelled correctly. Moreover, when looking for electromagnetic signatures of accretion (e.g., Tremou et al. 2018), it is crucial to have a good estimate for the wander radius of an IMBH within its host star cluster. An underestimate may lead offcentre radio sources, which could be potential IMBHs, to be excluded.
Finally, when an IMBH claim is made based on radial velocity signatures, as in the case of the IMBH in the Leo I dwarf spheroidal (BustamanteRosell et al. 2021), velocity dispersion anisotropy could be an important source of confusion (Zocchi et al. 2016), as strongly radially anisotropic systems could be compatible with a massive central object as well as with radially biased initial conditions.
Estimating the probability of IMBH expulsion and the wander radius, and tracing the evolution of anisotropy in the presence of an IMBH are three applications for which the recently introduced MPCDSS code (Di Cintio et al. 2021, 2022) becomes competitive in terms of realism with direct Nbody codes and other approximate approaches, as we argue in the following.
Direct Nbody simulations of stellar systems are often perceived as more realistic (e.g., see Takahashi & Portegies Zwart 2000; Baumgardt 2001; Hurley et al. 2005; Baumgardt et al. 2008; Bortolas et al. 2016; Wang 2020) than other approaches, such as Monte Carlo methods for instance (e.g., see Freitag & Benz 2001; Giersz 2006; Hypki & Giersz 2013; Giersz et al. 2013; Sollima & Mastrobuono Battisti 2014; Vasiliev 2015; Sollima & Ferraro 2019; Aros et al. 2020), especially when collisional dynamics is involved (see the discussions in Kim et al. 2008; Heggie 2011; Kamlah et al. 2022).
However, star clusters that are both in a collisional dynamic regime and contain N > 10^{6} particles are a common occurrence, even in our Galaxy (see the discussion in Di Cintio et al. 2021). Direct Nbody models of these star clusters cannot be simulated using a onetoone star to stellarparticle ratio, due to the computational constraints of the method. This greatly reduces the faithfulness of direct Nbody simulations.
When dealing with IMBH hosting systems, this limitation can be recast in terms of two dimensionless ratios that are important for determining the dynamical evolution of the system: the ratio of IMBH mass to the average stellar mass in the system,
and the ratio of IMBH mass to the total mass in the system
We can better appreciate the role of these two ratios by way of a very simplified example, where the star cluster has a typical size R, is virialized, and equipartition of kinetic energies holds between the IMBH and the surrounding stars with velocity dispersion σ. With these assumptions, the radius of the sphere of influence of the IMBH (i.e. the radius below which the BH potential Φ_{IMBH} = −GM_{IMBH}/r dominates over the contribution of the stellar component, e.g., see Peebles 1972; see also Merritt 2004) is
while the socalled wander radius (the typical distance at which the IMBH is found from the host centre of mass, see Bahcall & Wolf 1976; Brockamp et al. 2011) works out to
In other words, the sphere of influence is the region within which the motion of a star is heavily affected by the presence of the IMBH, while the wander radius is the typical distance at which we expect to find the IMBH from the bottom of the cluster potential well.
Seeing how r_{inf} and r_{wan} describe two different but equally important aspects of the IMBHhost interaction and how they depend on the two ratios introduced above, we must conclude that a simulation must match both α and μ of a real star cluster to be considered a realistic model thereof. However, the number of particles in a simulation is by definition
so a constraint on the number of particles N that can be simulated via direct Nbody simulations translates to an inaccessible region in the (μ, α) plane.
In Fig. 1 we show observational claims of IMBH detections in Milky Way globular clusters (GCs) in the (μ, α) plane. IMBH masses are based on Table 1 of Mezcua (2017), which includes maximum masses in case of a negative claim. The average stellar mass in each system is assumed to be 0.5 M_{⊙} and the total GC mass is taken from Baumgardt & Hilker (2018). As argued above, a correct modelling of the dynamical effects of an IMBH on its host GC must match both mass ratios. Figure 1 shows that this is simply not possible for about half of the claims reported by Mezcua (2017) unless simulations are run with at least a million stellar particles. This is at the limit of the current direct Nbody simulation state of the art.
Fig. 1. Observational claims of IMBH detection in the (μ = M_{IMBH}/⟨m⟩, α = M_{IMBH}/M) plane, shown as black diamonds. Any technique (most notably the direct Nbody one) that cannot simulate more than a given number of particles N must trade α for μ at any given N, being unable to access the top left region, which is shown shaded in orange for a typical value of N for direct Nbody codes currently available, N = 10^{6}. As a reference, N = 10^{5} is also shown as a dashed line. 
An example may help clarify the meaning of Fig. 1. If a simulation contains 10^{5} particles, and if the typical mean stellar mass for an old star cluster is 0.5 M_{⊙}, we can consider the simulation to model a 5 × 10^{4} M_{⊙} star cluster with a onetoone startoparticle ratio. If we wish to study the effect of an IMBH whose size is 0.5% of the total mass of this cluster, then the IMBH must have a mass of 2.5 × 10^{2} M_{⊙}, which is a rather small IMBH in comparison to the stars with which it interacts. Alternatively, we can imagine that the stellar particles in our simulation do not track actual stars onetoone, but then any claims that direct Nbody models are more realistic than other approximate methods would become untenable, and it is harder to justify applying stellar evolution to each stellar particle as if it were an individual star.
In this article we thus leverage our recently introduced MPCDSS code, which treats two and multiplebody collisions in an approximate fashion, to simulate largeN systems with a onetoone particletostar ratio while correctly modelling the IMBHhost interaction in the sense discussed above. In Di Cintio et al. (2022), we observed that while the presence of a mass spectrum generally speeds up the core collapse of a given model with respect to the parent system with the same total mass and N but with equal mass particles, a central IMBH typically induces a shallower core collapse. It has yet to be determine how the velocity dispersion and density profiles are affected by the presence of a central IMBH, and more relevantly how the anisotropy profile even evolves. Here we performed additional numerical experiments with a broader range of IMBH masses and different initial anisotropy profiles for equilibrium models with equal masses and the Salpeter (1955) initial mass function (IMF).
The paper is structured as follows: In Sect. 2 we discuss the simulation setup and the generation of the initial conditions, in addition to introducing the structure of MPCDSS. In Sect. 3 we present the results on the evolution of dynamical models of star clusters with a central IMBH. Finally, we provide a summary in Sect. refsec4.
2. Simulations
2.1. Initial conditions
In this work we have run a set of hybrid numerical simulations with our MPCDSS code discussed in DC2020, which combine a standard particlemesh approach for the stellar potential with a multiparticle collision scheme for the collisions to a direct Nbody code for the dynamics of the BH(s). The direct Nbody treatment of the BH is a negligible overhead with respect to pure MPCDSS since the BH is only one body.
We performed simulations with 3 × 10^{3} ≤ N ≤ 10^{6}, and adopted the usual Plummer (1911) profile
as an initial condition, with total mass M and scale radius r_{s} set to unity. Particle masses m_{i} are either equal to M/N or extracted from a Salpeter (1955) powerlaw mass function with the exponent −2.35 truncated such that the minimumtomaximummass ratio ℛ = m_{min}/m_{max} equals 10^{−3}.
In the runs without a central IMBH, we extracted the initial particles’ velocities using the rejection method on the numerically recovered phasespace distribution function with OsipkovMerritt (hereafter OM, Osipkov 1979; Merritt 1985) radial anisotropy defined by
In the equation above, , with ℰ and J being the particle’s energy and angular momentum per unit mass^{1}, respectively, Φ being the gravitational potential of the model, r_{a} being the anisotropy radius, and ρ_{a} being the augmented density, defined by
The anisotropy radius r_{a} is the control parameter associated with the extent of velocity anisotropy of the model, so that, for a given density profile, the velocitydispersion tensor is nearly isotropic inside r_{a}, and more and more radially anisotropic for increasing r. In other words, small values of r_{a} correspond to more radially anisotropic systems, and thus to larger values of the anisotropy parameter ξ (see e.g., Binney & Tremaine 2008) defined by
where K_{r} and K_{t} = K_{θ} + K_{ϕ} are the radial and tangential components of the kinetic energy tensor that read as
where, and are the radial and tangential phasespaceaveraged square velocity components, respectively.
For the Plummer density distribution (6), f(Q) is given explicitly in terms of elementary functions (e.g., see Dejonghe 1987; Breen et al. 2017) as
where is the (scalar) central velocity dispersion.
In the simulations featuring a central IMBH, again we extracted the particles’ positions from the density distribution (6). The correspondent velocities were sampled from the standard phasespace distribution given by Eq. (11) if their radial position r is larger than the influence radius r_{inf}, while instead for r ≤ r_{inf} the velocities were generated sampling the isotropic distribution function f(ℰ) for a homogeneous and noninteracting ‘atmosphere’ of density and radius r_{inf} in equilibrium in Φ_{IMBH}, which reads as
For the specific choice for a Plummer density profile, the influence radius is given as a function of the Plummer’s scale radius r_{s} and the IMBH mass in units of the mass ratio α as
We note that in previous work (e.g., see Chatterjee et al. 2002a,b), the initial conditions for the stellar component have always been set up by sampling the phasespace distribution for a Plummer model without the BH and later renormalized so that the resulting systems’ stars+BH system was virialized. We note also that, in principle, a system with a cored density profile (such as the Plummer used here) cannot have a consistent equilibrium phasespace distribution when embedded in an external potential associated with a singular density profile (such as that of the central BH, see Ciotti 1996).
In line with all of these works, we simulated the IMBH by adding a particle, initially sitting at rest in the centre of the system with mass 10^{−4}M ≤ M_{IMBH} ≤ 10^{−2}M. For the range of simulation particles 3 × 10^{3} ≤ N ≤ 10^{6}, such a choice corresponds to a range in mass ratios 3 ≤ μ ≤ 3 × 10^{3}. In Table 1 we summarize the parameter of the simulations discussed in the next sections.
Summary of the initial conditions.
2.2. Numerical scheme
Following Di Cintio et al. (2021, 2022), we evolved all sets of isolated initial conditions up to 2 × 10^{4} dynamical times , so that in all cases the systems reach core collapse and are evolved further after it for at least another 10^{3}t_{dyn}. We employed our recent implementation of MPCDSS where the gravitational potential and force were computed by the standard particleincell scheme on a fixed spherical grid of N_{g} = N_{r} × N_{ϑ} × N_{φ} mesh points (e.g., see Londrillo & Messina 1990).
In the simulations presented here, we have used N_{r} = 1024, N_{ϑ} = 16, and N_{φ} = 16 with logarithmically spaced radial bins and averaged the potential along the azimuthal and polar coordinates in order to enforce the spherical symmetry throughout the simulation.
The multiparticle collisions (see Di Cintio et al. 2017, 2021, 2022, 2023 for details) were performed on a different mesh with N_{g} = 32 × 16 × 16 extended only up to r_{cut} = 100r_{s} and conditioned with a standard rejection step to the local (i.e. cell dependent) collision probability p_{i} given by
where Δt is the simulation time step, ν_{c} is the collision frequency, β is a dimensionless constant of the order of twice the number of the simulation cells, and Erf(x) is the standard error function.
In Eq. (14), β is a dimensionless constant of the order of the total number of cells in the system and the collision frequency is defined as usual as
where n_{i} the local stellar number density, and σ_{i} the average particle mass and the putative velocity dispersion in the cell, respectively, and the Coulomb logarithm log Λ is fixed to 10.
In all simulations presented here, we used the same normalization such that G = M = r_{s} = t_{dyn} = v_{s} = 1. Hereafter (except where otherwise stated), all distances and velocities are given in units of the Plummer scale radius r_{s} and scale velocity v_{s} ≡ r_{s}/t_{dyn}. In our simulations, for such a choice for our units, we adopted a constant times step Δt and used a second order leap frog scheme to propagate the particle’s equations of motion. The specific (fixed) value of the time step in units of t_{dyn} depends on the number of particles N and their mean closest approach distance (e.g., see Dehnen & Read 2011) and ranges from 3 × 10^{−3} for N = 3 × 10^{3} to 1.25 × 10^{−2} for N = 10^{6}.
In the runs including the central IMBH, its interaction with the stars was evaluated directly, that is to say the IMBH did not take part in the MPC step nor in the evaluation of the mean field potential. In order to keep the same rather large Δt of the simulation, the potential exerted by the IMBH was regularized as
where we used ϵ = 10^{−4} in units of r_{s} so that for the IMBH masstocluster mass ratio 10^{−3}, the softening length was always of the order of one tenth of the influence radius of the IMBH r_{inf}. With such choices for simulation parameters, on average, the MPC simulations on a single core are a factor ∼10 faster than direct Nbody simulations for N of the order of 10^{4}, and they remain faster down to a factor ∼2 for N = 10^{6}.
3. Results
3.1. Evolution of density and central velocity dispersion
As an indicator of the evolution of the concentration of a given system, we followed the evolution of the Lagrangian radii containing a given fraction of the system’s mass and within such radii we also computed the mean velocity dispersion. In Fig. 2 (upper panels), we show the evolution of the Lagrangian radius r_{5%} enclosing 5% of a N = 3 × 10^{5} initially isotropic (i.e. ξ_{0} = 1) Plummer model with Salpeter IMF (left panels) and equal masses (right panels) and different values for the IMBH mass ratio α. As expected, models with a mass spectrum contract on shorter timescales (at least when the specific value for M_{IMBH} is low in units of M), with respect to their counterparts with equal mass particles. For large values of α, r_{5%} grows rapidly without showing signs of an earlier contraction. This implies that the presence of a massive IMBH should be associated with an inflated core (see also Fig. 9 in Di Cintio et al. 2022).
Fig. 2. Evolution of the Lagrangian radii r_{5%} enclosing 5% of the system’s mass (top panels) and central velocity dispersion σ_{c} evaluated for particles inside r_{5%} (bottom panels) for models with a Salpeter mass function (left) and equal masses (right), and α = 10^{−4}, 3 × 10^{−4}, 10^{−3}, 3 × 10^{−3}, and 10^{−2}. In all cases the systems were initially isotropic and N = 10^{5} so that μ = α × 10^{5}. 
The evolution of the average central velocity dispersion σ_{c} evaluated within r_{5%} is shown in the lower panels of Fig. 2 (see also the last column in Table 2). We observe that for the models with a mass spectrum, σ_{c} steadily decreases for the cases with α > 10^{−3}, while it grows reaching its maximum (surprisingly independent of α) at around t_{cc} and then decreases for the systems hosting a lowermass IMBH. In models in which all stars have the same mass, the behaviour of σ_{c} is the same for α > 10^{−3}, while it appears somewhat more complex at lower αs. Remarkably, σ_{c} reaches different maximum values before starting to decrease for different values of the mass ratio α. In other words, one can conclude that the presence of a massive IMBH in a star cluster, long after its core collapse, should induce a colder and larger core with respect to a star cluster in the same mass range but without a central IMBH.
Summary of the simulation properties.
3.2. Radial anisotropy profiles
Before exploring the effects of a central IMBH on the orbital anisotropy of a given model, we studied the evolution of the anisotropy profiles for systems without an IMBH well beyond core collapse (and mass segregation). For OM systems characterized by different values of N and ξ, we have evaluated the radial anisotropy profile (see e.g., Binney & Tremaine 2008)
at different times. We find that, surprisingly, for all values of N considered here between 3 × 10^{3} and 10^{6}, all isotropic models (ξ_{0} = 1, β(r, 0) = 0 everywhere) with a (Salpeter) mass spectrum have already evolved right before core collapse (typically at around ∼40t_{dyn}) in an ‘isotropic core’ for r < 3r_{s} surrounded by an increasingly anisotropic halo of weakly bound particles kicked out during the process of mass segregation. At later times, the profile of β remains relatively unchanged, as shown for the N = 10^{6} case in the left panel of Fig. 3 showing β at t = 100, 2000, and 10 000t_{dyn} (solid lines) and at the time of core collapse t_{cc} (thick dotteddashed line).
Fig. 3. Evolution of the radial anisotropy profile β(r) in Plummer models without a central IMBH, with a Salpeter IMF, N = 10^{6}, and from left to right ξ_{0} = 1 (isotropic), 1.5 (limit for stability), and 2.5 (critical, for consistency). The thin dashed lines mark the initial (analytical) anisotropy profile, while the dotdashed lines mark the profile β at the indicated time of core collapse t_{cc}. 
The systems starting with initial conditions sampled from OM models with a larger degree of anisotropy (cf. middle and right panels of Fig. 3) remarkably become less and less radially anisotropic, with respect to their initial state, marked in the figure by the thin dashed lines^{2}. We verified that such behaviour holds true even for other powerlaw mass spectra (not shown here) proportional to m^{−0.6}, m^{−1}, and m^{−3}. In these cases, the profile of β at t_{cc} for systems with low and large values for the mass function slope α are qualitatively very similar for both highly anisotropic Plummer initial conditions – which are closer to the critical value of the anisotropy indicator – for consistency (i.e. ξ_{0} = 2.5), as well as for moderately anisotropic initial conditions (i.e. ξ_{0} = 1.5); both cases show almost isotropic ‘cores’ up to r ≈ 10.
In practice, independently of the specific values of ξ_{0} > 1, the anisotropy radius of the model increases with time as the system undergoes core collapse and reexpands. Vice versa, initially isotropic star clusters do ‘anisotropize’ during core collapse and mass segregation, though their anisotropy radius also increases for t ≫ t_{cc}. Isotropic equal masses’ models (not shown here), with substantially longer core collapse timescales in the absence of a mass spectrum (see Col. 2 in Table 2), remain substantially isotropic everywhere (i.e. with final anisotropy radii usually at about 5r_{50%}), while OManisotropic equal masses’ models also experience an increase in anisotropy up to core collapse as their multimass counterparts. We observe that, in general, the models have a longer core collapse timescale at a fixed mass for increasing values of the initial anisotropy parameter ξ_{0}, independently of the specific mass spectrum.
Adding a central BH, for all of the initial conditions discussed here, has the effect of systematically reducing t_{cc} of a factor between 2 and 3.5. As observed for models without a central BH, isotropic initial conditions tend to evolve towards more anisotropic states also for the cases with a BH (see left panels in Fig. 4) with or without a mass spectrum. In the latter case, the strong kicks exerted by the BH on eccentric orbits play the role of mass segregation in populating the outer radii of low angular momentum stars. OM models with a central BH again evolve towards less anisotropic states with equal mass systems with significantly flatter β profiles at late times. Of course, the interplay between the evolution of the anisotropy profiles and that of the IMBH should be, in principle, studied in models starting from initial conditions where the protocluster is far from being virialized with or without a significantly massive seed for the IMBH, such as those produced in Torniamenti et al. (2022).
Fig. 4. Evolution of the radial anisotropy profile β(r) in Plummer models with a central IMBH of mass M_{IMBH} = 10^{−3}M, with a Salpeter IMF (top panel row) or equal masses (bottom panel row), N = 10^{6}, and from left to right ξ_{0} = 1 (isotropic), 1.5 (limit for stability), and 2.5 (critical for consistency). The thin dashed lines mark the initial (analytical) anisotropy profile. The empirical lines become somewhat noisy to the left because of the increasingly low number of particles available for calculating β(r) at small radii. 
3.3. Wander radius
For the IMBH hosted in the models discussed in the previous sections, we have evaluated the probability density functions (PDF) of the radial position with respect to the geometric centre of the stellar distribution f(r) and velocity f(v). In Fig. 5, we show said distributions for initially isotropic models with N = 10^{6} and 10^{−4} ≤ α ≤ 10^{−2} and either singlemass or Salpeter mass spectra, respectively. As an indicator of the extension of the BH wander radius r_{wan}, we extracted the radius corresponding to the peak of f(r). We observe that (see Fig. 6) for fixed N while the peak of the velocity distribution f(v) moves at lower velocities for increasing α (for both equal mass and Salpeter systems), the peak of the radial position distribution f(r) (i.e. the putative wander radius) is somewhat independent of α for equal mass models, while it moves to smaller values of r for increasing α in systems with a mass spectrum.
Fig. 5. Radial and velocity distributions for the IMBHs. Upper panels: distributions of the radial coordinate (left) and velocity (right) attained by IMBHs hosted in a Plummer model with N = 10^{5} equal mass particles with an initially isotropic velocity distribution and the same combinations of α and μ as in Fig. 2. Lower panels: same as above but for models with Salpeter MF. 
Fig. 6. Most probable velocity and putative wander radius of the IMBH, respectively, shown as a function of the mass ratio α for single mass models (empty circles) and models with a Salpeter mass function (filled circles). In all cases, N = 10^{5}. 
When fixing α while decreasing μ (i.e. we changed N so that the BH masstomean stellar mass ⟨m⟩ varied), both f(r) and f(v) become broader and peak at larger r and v, respectively (see Fig. 7, main panels). Remarkably, we recovered (at least for μ ≲ 200) the predicted r_{wan} ∝ μ^{−1/2} trend, while we observe a clear v_{IMBH} ∝ μ^{−1/3} behaviour for the typical velocity of the IMBH (dashed and dotted lines in the right panel of Fig. 7, respectively). Not surprisingly, for fixed μ we observe larger values for the wander radius r_{wan} in systems with a smaller α, as in those cases the larger cluster mass forms a deeper potential well. This is exemplified in Fig. 8 where r_{wan} (top panel) and v_{IMBH} (bottom panel) are plotted and colour coded against α and μ for the equal m cases.
Fig. 7. Distributions of the radial coordinate (left panel) and velocity (middle panel) of the IMBHs embedded in a Plummer model with 3 × 10^{3} ≤ N ≤ 10^{6} equal mass particles and an initially isotropic velocity distribution. The top right and top left panels show the peaks of the velocity and radial distributions as function of the mass ratio μ, respectively. Furthermore, shows a markedly μ^{−1/3} trend (dotted line), while r_{wan} has a trend compatible with μ^{−1/2} for μ < 3 × 10^{2} (dash ed line). In all cases, α = 10^{−3}. 
Fig. 8. Wander radius (scaled to 01 over the simulation sample) as a function of μ and α (top panel) and wander velocity as a function of μ and α (bottom panel). The diagonal lines correspond to a constant number of particles: 10^{4} (solid), 10^{5} (dashed), 10^{6} (dotted), and 10^{7} (dotdashed). 
We observe that, in general, for mass ratios μ larger than ten, all such trends are weakly affected by the mass spectrum or the specific anisotropy profile of the model at hand. However, we notice that for increasing initial values of ξ, the distribution of the radial coordinate of the IMBH shows systematically fatter tails, corresponding to a decreasing (negative) kurtosis κ. As an example, in Fig. 9 we show f(r) for μ = 10^{3}, α = 10^{−3}, and ξ_{0} = 1, 1.5, and 2.5, as well as Salpeter (left panel) and equal mass (right panel) models. This implies that IMBHs in models with markedly anisotropic initial conditions might have a nonnegligible probability of being displaced by a few scale radii from the geometric centre of the star cluster. We note that Chatterjee et al. (2002a,b), by means of direct Nbody simulations, and FokkerPlanck calculations in a static cluster potential Φ estimated a limit r_{wan} of the order of 0.1r_{s}, where r_{s} is some scale length roughly equal to the halfmass radius of the model at hand. Such a value is typically assumed as the radial distance^{3} within which to look for IMBH candidates in GCs in many observational studies. It is important to note that, on the one hand, such Nbody runs had the natural limits of the relatively small number of (equal mass) particles N and their overall large computational cost. On the other hand, the FokkerPlanck models used, somewhat arbitrarily, a Gaussian force fluctuation distribution f(δF). Therefore, in both cases, the rare but strong encounters were systematically neglected.
Fig. 9. Distribution of the IMBH radial coordinate for Plummer models with N = 10^{6}, M_{IMBH} = 10^{−3}, and ξ = 1, 1.5, and 2.5, for the model with a Salpeter mass function (left panel) and equal masses (right panel). 
Di Cintio et al. (2020) studied the dynamics of massive BH in galactic cores using a model based on the integration of stochastic (i.e. Langevin) equations (see also Pasquato & Di Cintio 2020) of the form
where η is the Chandrasekhar dynamical friction coefficient and δF is a fluctuating force (per unit mass). They show that the position distribution of the BH extracted from short time Nbody simulations is qualitatively ‘intermediate’ between those obtained in longer Langevin simulations with force fluctuations that are sampled from a Gaussian and a Holtsmark (1919) distribution for f(δF) (see their Fig. 1). We stress the fact that the Holtsmark distribution is the correct one that models the force fluctuations in a system of particles interacting with a 1/r^{2} force law (Chandrasekhar & von Neumann 1942, 1943). We also remark that, in the MPC simulations discussed in the present work the interactions between the IMBH and the stars are evaluated with a direct sum scheme (as in Chatterjee et al. 2002a), but the evolution time, being of the order of several thousands of crossing times, is much larger, thus allowing for strong encounters (typically corresponding to strong force fluctuations described by the heavy F^{−5/2} tails of the Holtsmark distribution) to have a nonnegligible role in the dynamics of the IMBH.
3.4. Escapers’ and compact objects’ retention fraction
In Di Cintio et al. (2021), we compare the timedependent fraction of escapers (i.e. particles reaching, with positive total energy, a truncation radius fixed at ∼20r_{50%}) in MPC and direct simulations with isotropic initial conditions with N of the order of 10^{4}, finding a rather good agreement for several choices for the mass spectrum. Here we evaluate the fraction of escapers for a broader range of N and different choices for ξ_{0} (cf. Col. 3 in Table 2).
In general, over a time span of 10^{4}t_{dyn}, models with equal masses tend to have a lower fraction of escapers than those with the same N and ξ_{0} with a mass spectrum; this is ascribed to the mass segregation process that pushes heavier stars to the inner regions of the cluster lowering their potential energy, at the expense of lighter stars that are pushed outside with increasing kinetic energies. This is also observable in Di Cintio et al. (2021; cf. Fig. 6 therein), where models characterized by heavier tailed mass functions (i.e. larger fractions of heavy particles at fixed ⟨m⟩) show a steeper time increase in the escapers’ fraction.
The presence of an IMBH, even if associated with a shallower core collapse, typically enhances particle evaporation via direct collisions with larger escapers’ fractions for increasing values of the mass ratio μ. For the models with or without an IMBH, the initial anisotropy profile has little influence on the fraction of escapers and no apparent trend is evident, as the latter might depend on μ, ξ, and α simultaneously (see Fig. 10 below).
Fig. 10. Fraction of escapers as a function of μ and α (top panel) as a function of μ and α (bottom panel). The diagonal lines correspond to a constant number of particles: 10^{4} (solid), 10^{5} (dashed), 10^{6} (dotted), and 10^{7} (dotdashed). 
For the models with a Salpeter mass function, we have also evaluated the massdependent escaper fraction, dividing the mass spectrum of the system into 50 logarithmically spaced mass bins. In Col. 4 of Table 2, we give the fraction of escapers in the largest bin (corresponding roughly to m > 25⟨m⟩). For star clusters with a mean stellar mass of about 0.5 M_{⊙}, these would correspond to m > 10 M_{⊙}, likely encompassing collapsed objects. Not surprisingly, bigger values of μ are associated with increasing fractions of heavy escapers. In the worst case, up to 35% of particles in the largest mass bin are ejected before 10^{4}t_{dyn}, which corresponds to a compact object retention fraction of about 65%. Again, the initial anisotropy profile does not have a significant effect on the retention fraction for fixed values of μ or α.
4. Discussion and conclusions
With multiparticle collision simulations with MPCDSS, we have investigated the dynamics of IMBHs in star clusters under different characteristics of the host (mass spectrum, orbital anisotropy) and the IMBH itself (mass ratio to the typical star and to the total host mass). Thanks to the linear complexity of MPCDSS with the number of particles, we have had the opportunity to explore a wider range in these mass ratios as discussed in the Introduction.
We have confirmed our preliminary results of Di Cintio et al. (2022) that the presence of a central BH with a mass of about 10^{−3} in units of the total cluster mass induces a faster but shallower core collapse. This remains true for other values of the mass rations μ and α defined in Sect. 1.
In practice, clusters hosting a central IMBH would be observed as ‘dynamically older’ than their counterparts with no BH and with a more diffuse and colder core. Moreover, for fixed mean stellar mass, all systems with the central BH have a significantly larger fraction of escapers (and a smaller retention fraction of heavier stars) than those with no BH. We have explored the effect of OsipkovMerritt initial anisotropy profiles, finding that long after the core collapse time was reached, independently of the initial value of ξ and the presence or lack thereof of an IMBH, the clusters show a anisotropy profile with β between two and ten initial scale radii r_{s}, or of about five final halfmass radii r_{50%}.
We have evaluated the PDF of the radial displacement of the IMBH f(r) (i.e. the distribution of its radial distance from the geometric centre of the star cluster) and defined its absolute maximum as an IMBH wander radius. If on the one hand, we recover the (M_{IMBH}/⟨m⟩)^{−1/2} trend (independently of the cluster mass and anisotropy profile), on the other hand, we observe that such a radius is seemingly less dependent on the α ratio, being typically of the order of 10^{−1}r_{50%}. A result whose importance cannot be overstated is that the distribution of the distance as well as the skyprojected distance from the centre attained by the IMBH in our simulations becomes distinctly leptokurtic for increasing values of the systems’ initial anisotropy. This corresponds to the presence of heavy tails, with the associated risk of underestimating the probability of low probability events. A possible astrophysical consequence could be the unduly exclusion of potential IMBH candidates when they happen to be too far away from the host systems’ centre based on our Gaussian and Brownian expectations. Tremou et al. (2018), for instance, exclude several radio sources from their analysis even though they are relatively near to the host star cluster centre, because they are further out than the estimated Brownian radius of an IMBH of the relevant mass.
We note that, by doing so we are assuming that the degree of anisotropy is independent of the specific particle mass. In principle, it would also be possible to generate initial conditions where different masses are associated with different degrees of radial anisotropy, e.g., see Gieles & Zocchi (2015), and also Webb et al. (2023).
Such a radius of about 0.1r_{s} is also consistent with the typical corestalling radius where dynamical friction and dynamical buoyancy compensate for each other (e.g., see Banik & van den Bosch 2021, 2022).
Acknowledgments
This material is based upon work supported by the “Fondazione Cassa di Risparmio di Firenze” under the project HIPERCRHEL for the use of high performance computing resources at the university of Firenze. P.F.D.C. is supported by the MIURPRIN2017 project Coarsegrained description for nonequilibrium systems and transport phenomena (CONEST) n.201798CZL. L.B. is financed by the “Fondazione Cassa di Risparmio di Firenze” under the project THE SWITCH. M. P. acknowledges financial support from the European Union’s Horizon 2020 research and innovation programme under the Marie SklodowskaCurie grant agreement No. 896248. We warmly acknowledge Evangelia Tremou for reading a draft of this manuscript.
References
 Aros, F. I., Sippel, A. C., MastrobuonoBattisti, A., et al. 2020, MNRAS, 499, 4646 [CrossRef] [Google Scholar]
 Bahcall, J. N., & Wolf, R. A. 1976, ApJ, 209, 214 [NASA ADS] [CrossRef] [Google Scholar]
 Ballone, A., Mapelli, M., & Pasquato, M. 2018, MNRAS, 480, 4684 [CrossRef] [Google Scholar]
 Banik, U., & van den Bosch, F. C. 2021, ApJ, 912, 43 [NASA ADS] [CrossRef] [Google Scholar]
 Banik, U., & van den Bosch, F. C. 2022, ApJ, 926, 215 [NASA ADS] [CrossRef] [Google Scholar]
 Baumgardt, H. 2001, MNRAS, 325, 1323 [NASA ADS] [CrossRef] [Google Scholar]
 Baumgardt, H., & Hilker, M. 2018, MNRAS, 478, 1520 [Google Scholar]
 Baumgardt, H., De Marchi, G., & Kroupa, P. 2008, ApJ, 685, 247 [NASA ADS] [CrossRef] [Google Scholar]
 Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition (Princeton University Press) [Google Scholar]
 Bortolas, E., Gualandris, A., Dotti, M., Spera, M., & Mapelli, M. 2016, MNRAS, 461, 1023 [Google Scholar]
 Breen, P. G., Varri, A. L., & Heggie, D. C. 2017, MNRAS, 471, 2778 [NASA ADS] [CrossRef] [Google Scholar]
 Brockamp, M., Baumgardt, H., & Kroupa, P. 2011, MNRAS, 418, 1308 [CrossRef] [Google Scholar]
 BustamanteRosell, M. J., Noyola, E., Gebhardt, K., et al. 2021, ApJ, 921, 107 [NASA ADS] [CrossRef] [Google Scholar]
 Carr, B. J., Bond, J. R., & Arnett, W. D. 1984, ApJ, 277, 445 [NASA ADS] [CrossRef] [Google Scholar]
 Chandrasekhar, S., & von Neumann, J. 1942, ApJ, 95, 489 [NASA ADS] [CrossRef] [Google Scholar]
 Chandrasekhar, S., & von Neumann, J. 1943, ApJ, 97, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Chatterjee, P., Hernquist, L., & Loeb, A. 2002a, Phys. Rev. Lett., 88, 121103 [NASA ADS] [CrossRef] [Google Scholar]
 Chatterjee, P., Hernquist, L., & Loeb, A. 2002b, ApJ, 572, 371 [NASA ADS] [CrossRef] [Google Scholar]
 Ciotti, L. 1996, ApJ, 471, 68 [NASA ADS] [CrossRef] [Google Scholar]
 Dehnen, W., & Read, J. I. 2011, Eur. Phys. J. Plus, 126, 55 [Google Scholar]
 Dejonghe, H. 1987, MNRAS, 224, 13 [NASA ADS] [CrossRef] [Google Scholar]
 Di Carlo, U. N., Mapelli, M., Pasquato, M., et al. 2021, MNRAS, 507, 5132 [NASA ADS] [CrossRef] [Google Scholar]
 Di Cintio, P., Livi, R., Lepri, S., & Ciraolo, G. 2017, Phys. Rev. E, 95, 043203 [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 Cintio, P., Pasquato, M., Kim, H., & Yoon, S.J. 2021, A&A, 649, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Di Cintio, P., Pasquato, M., SimonPetit, A., & Yoon, S.J. 2022, A&A, 659, A19 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Di Cintio, P., Pasquato, M., Barbieri, L., et al. 2023, IAU Symp., 362, 134 [NASA ADS] [Google Scholar]
 Ebisuzaki, T., Makino, J., Tsuru, T. G., et al. 2001, ApJ, 562, L19 [NASA ADS] [CrossRef] [Google Scholar]
 Freitag, M., & Benz, W. 2001, A&A, 375, 711 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gieles, M., & Zocchi, A. 2015, MNRAS, 454, 576 [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]
 Greene, J. E., Strader, J., & Ho, L. C. 2020, ARA&A, 58, 257 [Google Scholar]
 Greif, T. H. 2015, Comput. Astrophys. Cosmol., 2, 3 [NASA ADS] [CrossRef] [Google Scholar]
 Heggie, D. C. 2011, Bull. Astron. Soc. India, 39, 69 [NASA ADS] [Google Scholar]
 HolleyBockelmann, K., Gültekin, K., Shoemaker, D., & Yunes, N. 2008, ApJ, 686, 829 [Google Scholar]
 Holtsmark, J. 1919, Ann. Phys., 363, 577 [CrossRef] [Google Scholar]
 Hurley, J. R., Pols, O. R., Aarseth, S. J., & Tout, C. A. 2005, MNRAS, 363, 293 [NASA ADS] [CrossRef] [Google Scholar]
 Hypki, A., & Giersz, M. 2013, MNRAS, 429, 1221 [NASA ADS] [CrossRef] [Google Scholar]
 Inayoshi, K., Visbal, E., & Haiman, Z. 2020, ARA&A, 58, 27 [NASA ADS] [CrossRef] [Google Scholar]
 Kamlah, A. W. H., Leveque, A., Spurzem, R., et al. 2022, MNRAS, 511, 4060 [NASA ADS] [CrossRef] [Google Scholar]
 Kim, E., Yoon, I., Lee, H. M., & Spurzem, R. 2008, MNRAS, 383, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Li, G.P. 2022, A&A, 666, A194 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lodato, G., & Natarajan, P. 2006, MNRAS, 371, 1813 [NASA ADS] [CrossRef] [Google Scholar]
 Loeb, A., & Rasio, F. A. 1994, ApJ, 432, 52 [Google Scholar]
 Londrillo, P., & Messina, A. 1990, MNRAS, 242, 595 [NASA ADS] [CrossRef] [Google Scholar]
 Merritt, D. 1985, AJ, 90, 1027 [Google Scholar]
 Merritt, D. 2004, in Coevolution of Black Holes and Galaxies, ed. L. C. Ho, 263 [Google Scholar]
 Mezcua, M. 2017, Int. J. Mod. Phys. D, 26, 1730021 [Google Scholar]
 Miller, M. C., & Hamilton, D. P. 2002, MNRAS, 330, 232 [Google Scholar]
 Osipkov, L. P. 1979, Soviet. Astron. Lett., 5, 42 [NASA ADS] [Google Scholar]
 Pacucci, F., & Loeb, A. 2022, MNRAS, 509, 1885 [Google Scholar]
 Pasquato, M., & Di Cintio, P. 2020, A&A, 640, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Peebles, P. J. E. 1972, ApJ, 178, 371 [NASA ADS] [CrossRef] [Google Scholar]
 Plummer, H. C. 1911, MNRAS, 71, 460 [Google Scholar]
 Portegies Zwart, S. F., Baumgardt, H., Hut, P., Makino, J., & McMillan, S. L. W. 2004, Nature, 428, 724 [Google Scholar]
 Rizzuto, F. P., Naab, T., Spurzem, R., et al. 2021, MNRAS, 501, 5257 [NASA ADS] [CrossRef] [Google Scholar]
 Rizzuto, F. P., Naab, T., Rantala, A., et al. 2023, MNRAS, 521, 2930 [NASA ADS] [CrossRef] [Google Scholar]
 Salpeter, E. E. 1955, ApJ, 121, 161 [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]
 Takahashi, K., & Portegies Zwart, S. F. 2000, ApJ, 535, 759 [NASA ADS] [CrossRef] [Google Scholar]
 Torniamenti, S., Pasquato, M., Di Cintio, P., et al. 2022, MNRAS, 510, 2097 [Google Scholar]
 Tremou, E., Strader, J., Chomiuk, L., et al. 2018, ApJ, 862, 16 [Google Scholar]
 Vasiliev, E. 2015, MNRAS, 446, 3150 [NASA ADS] [CrossRef] [Google Scholar]
 Volonteri, M., Lodato, G., & Natarajan, P. 2008, MNRAS, 383, 1079 [Google Scholar]
 Wang, L. 2020, MNRAS, 491, 2413 [NASA ADS] [Google Scholar]
 Weatherford, N. C., Kıroğlu, F., Fragione, G., et al. 2023, ApJ, 946, 104 [NASA ADS] [CrossRef] [Google Scholar]
 Webb, J. J., Hunt, J. A. S., & Bovy, J. 2023, MNRAS, 521, 3898 [NASA ADS] [CrossRef] [Google Scholar]
 Yoshida, N., Omukai, K., Hernquist, L., & Abel, T. 2006, ApJ, 652, 6 [NASA ADS] [CrossRef] [Google Scholar]
 Zocchi, A., Gieles, M., & HénaultBrunet, V. 2016, in Star Clusters and Black Holes in Galaxies across Cosmic Time, eds. Y. Meiron, S. Li, F. K. Liu, & R. Spurzem, 312, 197 [NASA ADS] [Google Scholar]
All Tables
All Figures
Fig. 1. Observational claims of IMBH detection in the (μ = M_{IMBH}/⟨m⟩, α = M_{IMBH}/M) plane, shown as black diamonds. Any technique (most notably the direct Nbody one) that cannot simulate more than a given number of particles N must trade α for μ at any given N, being unable to access the top left region, which is shown shaded in orange for a typical value of N for direct Nbody codes currently available, N = 10^{6}. As a reference, N = 10^{5} is also shown as a dashed line. 

In the text 
Fig. 2. Evolution of the Lagrangian radii r_{5%} enclosing 5% of the system’s mass (top panels) and central velocity dispersion σ_{c} evaluated for particles inside r_{5%} (bottom panels) for models with a Salpeter mass function (left) and equal masses (right), and α = 10^{−4}, 3 × 10^{−4}, 10^{−3}, 3 × 10^{−3}, and 10^{−2}. In all cases the systems were initially isotropic and N = 10^{5} so that μ = α × 10^{5}. 

In the text 
Fig. 3. Evolution of the radial anisotropy profile β(r) in Plummer models without a central IMBH, with a Salpeter IMF, N = 10^{6}, and from left to right ξ_{0} = 1 (isotropic), 1.5 (limit for stability), and 2.5 (critical, for consistency). The thin dashed lines mark the initial (analytical) anisotropy profile, while the dotdashed lines mark the profile β at the indicated time of core collapse t_{cc}. 

In the text 
Fig. 4. Evolution of the radial anisotropy profile β(r) in Plummer models with a central IMBH of mass M_{IMBH} = 10^{−3}M, with a Salpeter IMF (top panel row) or equal masses (bottom panel row), N = 10^{6}, and from left to right ξ_{0} = 1 (isotropic), 1.5 (limit for stability), and 2.5 (critical for consistency). The thin dashed lines mark the initial (analytical) anisotropy profile. The empirical lines become somewhat noisy to the left because of the increasingly low number of particles available for calculating β(r) at small radii. 

In the text 
Fig. 5. Radial and velocity distributions for the IMBHs. Upper panels: distributions of the radial coordinate (left) and velocity (right) attained by IMBHs hosted in a Plummer model with N = 10^{5} equal mass particles with an initially isotropic velocity distribution and the same combinations of α and μ as in Fig. 2. Lower panels: same as above but for models with Salpeter MF. 

In the text 
Fig. 6. Most probable velocity and putative wander radius of the IMBH, respectively, shown as a function of the mass ratio α for single mass models (empty circles) and models with a Salpeter mass function (filled circles). In all cases, N = 10^{5}. 

In the text 
Fig. 7. Distributions of the radial coordinate (left panel) and velocity (middle panel) of the IMBHs embedded in a Plummer model with 3 × 10^{3} ≤ N ≤ 10^{6} equal mass particles and an initially isotropic velocity distribution. The top right and top left panels show the peaks of the velocity and radial distributions as function of the mass ratio μ, respectively. Furthermore, shows a markedly μ^{−1/3} trend (dotted line), while r_{wan} has a trend compatible with μ^{−1/2} for μ < 3 × 10^{2} (dash ed line). In all cases, α = 10^{−3}. 

In the text 
Fig. 8. Wander radius (scaled to 01 over the simulation sample) as a function of μ and α (top panel) and wander velocity as a function of μ and α (bottom panel). The diagonal lines correspond to a constant number of particles: 10^{4} (solid), 10^{5} (dashed), 10^{6} (dotted), and 10^{7} (dotdashed). 

In the text 
Fig. 9. Distribution of the IMBH radial coordinate for Plummer models with N = 10^{6}, M_{IMBH} = 10^{−3}, and ξ = 1, 1.5, and 2.5, for the model with a Salpeter mass function (left panel) and equal masses (right panel). 

In the text 
Fig. 10. Fraction of escapers as a function of μ and α (top panel) as a function of μ and α (bottom panel). The diagonal lines correspond to a constant number of particles: 10^{4} (solid), 10^{5} (dashed), 10^{6} (dotted), and 10^{7} (dotdashed). 

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.