Open Access
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/0004-6361/202346124
Published online 24 April 2023

© The Authors 2023

Licence Creative CommonsOpen 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., Holley-Bockelmann 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 off-centre 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 (Bustamante-Rosell 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 N-body codes and other approximate approaches, as we argue in the following.

Direct N-body 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 > 106 particles are a common occurrence, even in our Galaxy (see the discussion in Di Cintio et al. 2021). Direct N-body models of these star clusters cannot be simulated using a one-to-one star to stellar-particle ratio, due to the computational constraints of the method. This greatly reduces the faithfulness of direct N-body 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,

(1)

and the ratio of IMBH mass to the total mass in the system

(2)

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 = −GMIMBH/r dominates over the contribution of the stellar component, e.g., see Peebles 1972; see also Merritt 2004) is

(3)

while the so-called 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

(4)

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 rinf and rwan describe two different but equally important aspects of the IMBH-host 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

(5)

so a constraint on the number of particles N that can be simulated via direct N-body 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 N-body simulation state of the art.

thumbnail Fig. 1.

Observational claims of IMBH detection in the (μ = MIMBH/⟨m⟩, α = MIMBH/M) plane, shown as black diamonds. Any technique (most notably the direct N-body 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 N-body codes currently available, N = 106. As a reference, N = 105 is also shown as a dashed line.

An example may help clarify the meaning of Fig. 1. If a simulation contains 105 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 × 104M star cluster with a one-to-one star-to-particle 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 × 102 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 one-to-one, but then any claims that direct N-body 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 multiple-body collisions in an approximate fashion, to simulate large-N systems with a one-to-one particle-to-star ratio while correctly modelling the IMBH-host 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 particle-mesh approach for the stellar potential with a multi-particle collision scheme for the collisions to a direct N-body code for the dynamics of the BH(s). The direct N-body 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 × 103 ≤ N ≤ 106, and adopted the usual Plummer (1911) profile

(6)

as an initial condition, with total mass M and scale radius rs set to unity. Particle masses mi are either equal to M/N or extracted from a Salpeter (1955) power-law mass function with the exponent −2.35 truncated such that the minimum-to-maximum-mass ratio ℛ = mmin/mmax equals 10−3.

In the runs without a central IMBH, we extracted the initial particles’ velocities using the rejection method on the numerically recovered phase-space distribution function with Osipkov-Merritt (hereafter OM, Osipkov 1979; Merritt 1985) radial anisotropy defined by

(7)

In the equation above, , with ℰ and J being the particle’s energy and angular momentum per unit mass1, respectively, Φ being the gravitational potential of the model, ra being the anisotropy radius, and ρa being the augmented density, defined by

(8)

The anisotropy radius ra is the control parameter associated with the extent of velocity anisotropy of the model, so that, for a given density profile, the velocity-dispersion tensor is nearly isotropic inside ra, and more and more radially anisotropic for increasing r. In other words, small values of ra correspond to more radially anisotropic systems, and thus to larger values of the anisotropy parameter ξ (see e.g., Binney & Tremaine 2008) defined by

(9)

where Kr and Kt = Kθ + Kϕ are the radial and tangential components of the kinetic energy tensor that read as

(10)

where, and are the radial and tangential phase-space-averaged 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

(11)

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 phase-space distribution given by Eq. (11) if their radial position r is larger than the influence radius rinf, while instead for r ≤ rinf the velocities were generated sampling the isotropic distribution function f(ℰ) for a homogeneous and non-interacting ‘atmosphere’ of density and radius rinf in equilibrium in ΦIMBH, which reads as

(12)

For the specific choice for a Plummer density profile, the influence radius is given as a function of the Plummer’s scale radius rs and the IMBH mass in units of the mass ratio α as

(13)

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 phase-space 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 phase-space 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−4M ≤ MIMBH ≤ 10−2M. For the range of simulation particles 3 × 103 ≤ N ≤ 106, such a choice corresponds to a range in mass ratios 3 ≤ μ ≤ 3 × 103. In Table 1 we summarize the parameter of the simulations discussed in the next sections.

Table 1.

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 × 104 dynamical times , so that in all cases the systems reach core collapse and are evolved further after it for at least another 103tdyn. We employed our recent implementation of MPCDSS where the gravitational potential and force were computed by the standard particle-in-cell scheme on a fixed spherical grid of Ng = Nr × Nϑ × Nφ mesh points (e.g., see Londrillo & Messina 1990).

In the simulations presented here, we have used Nr = 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 multi-particle collisions (see Di Cintio et al. 2017, 2021, 2022, 2023 for details) were performed on a different mesh with Ng = 32 × 16 × 16 extended only up to rcut = 100rs and conditioned with a standard rejection step to the local (i.e. cell dependent) collision probability pi given by

(14)

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

(15)

where ni 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 = rs = tdyn = vs = 1. Hereafter (except where otherwise stated), all distances and velocities are given in units of the Plummer scale radius rs and scale velocity vs ≡ rs/tdyn. 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 tdyn 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 × 103 to 1.25 × 10−2 for N = 106.

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

(16)

where we used ϵ = 10−4 in units of rs so that for the IMBH mass-to-cluster mass ratio 10−3, the softening length was always of the order of one tenth of the influence radius of the IMBH rinf. With such choices for simulation parameters, on average, the MPC simulations on a single core are a factor ∼10 faster than direct N-body simulations for N of the order of 104, and they remain faster down to a factor ∼2 for N = 106.

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 r5% enclosing 5% of a N = 3 × 105 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 MIMBH is low in units of M), with respect to their counterparts with equal mass particles. For large values of α, r5% 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).

thumbnail Fig. 2.

Evolution of the Lagrangian radii r5% enclosing 5% of the system’s mass (top panels) and central velocity dispersion σc evaluated for particles inside r5% (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 = 105 so that μ = α × 105.

The evolution of the average central velocity dispersion σc evaluated within r5% 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 tcc and then decreases for the systems hosting a lower-mass 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.

Table 2.

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)

(17)

at different times. We find that, surprisingly, for all values of N considered here between 3 × 103 and 106, all isotropic models (ξ0 = 1, β(r, 0) = 0 everywhere) with a (Salpeter) mass spectrum have already evolved right before core collapse (typically at around ∼40tdyn) in an ‘isotropic core’ for r < 3rs 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 = 106 case in the left panel of Fig. 3 showing β at t = 100, 2000, and 10 000tdyn (solid lines) and at the time of core collapse tcc (thick dotted-dashed line).

thumbnail Fig. 3.

Evolution of the radial anisotropy profile β(r) in Plummer models without a central IMBH, with a Salpeter IMF, N = 106, 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 dot-dashed lines mark the profile β at the indicated time of core collapse tcc.

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 lines2. We verified that such behaviour holds true even for other power-law mass spectra (not shown here) proportional to m−0.6, m−1, and m−3. In these cases, the profile of β at tcc 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 re-expands. Vice versa, initially isotropic star clusters do ‘anisotropize’ during core collapse and mass segregation, though their anisotropy radius also increases for t ≫ tcc. 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 5r50%), while OM-anisotropic equal masses’ models also experience an increase in anisotropy up to core collapse as their multi-mass 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 tcc 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 proto-cluster is far from being virialized with or without a significantly massive seed for the IMBH, such as those produced in Torniamenti et al. (2022).

thumbnail Fig. 4.

Evolution of the radial anisotropy profile β(r) in Plummer models with a central IMBH of mass MIMBH = 10−3M, with a Salpeter IMF (top panel row) or equal masses (bottom panel row), N = 106, 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 = 106 and 10−4 ≤ α ≤ 10−2 and either single-mass or Salpeter mass spectra, respectively. As an indicator of the extension of the BH wander radius rwan, 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.

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

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

When fixing α while decreasing μ (i.e. we changed N so that the BH mass-to-mean 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 rwan ∝ μ−1/2 trend, while we observe a clear vIMBH ∝ μ−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 rwan 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 rwan (top panel) and vIMBH (bottom panel) are plotted and colour coded against α and μ for the equal m cases.

thumbnail Fig. 7.

Distributions of the radial coordinate (left panel) and velocity (middle panel) of the IMBHs embedded in a Plummer model with 3 × 103 ≤ N ≤ 106 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 rwan has a trend compatible with μ−1/2 for μ < 3 × 102 (dash ed line). In all cases, α = 10−3.

thumbnail Fig. 8.

Wander radius (scaled to 0-1 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: 104 (solid), 105 (dashed), 106 (dotted), and 107 (dot-dashed).

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 μ = 103, α = 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 non-negligible 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 N-body simulations, and Fokker-Planck calculations in a static cluster potential Φ estimated a limit rwan of the order of 0.1rs, where rs is some scale length roughly equal to the half-mass radius of the model at hand. Such a value is typically assumed as the radial distance3 within which to look for IMBH candidates in GCs in many observational studies. It is important to note that, on the one hand, such N-body 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 Fokker-Planck models used, somewhat arbitrarily, a Gaussian force fluctuation distribution f(δF). Therefore, in both cases, the rare but strong encounters were systematically neglected.

thumbnail Fig. 9.

Distribution of the IMBH radial coordinate for Plummer models with N = 106, MIMBH = 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

(18)

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 N-body 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/r2 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 non-negligible role in the dynamics of the IMBH.

3.4. Escapers’ and compact objects’ retention fraction

In Di Cintio et al. (2021), we compare the time-dependent fraction of escapers (i.e. particles reaching, with positive total energy, a truncation radius fixed at ∼20r50%) in MPC and direct simulations with isotropic initial conditions with N of the order of 104, 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 104tdyn, 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).

thumbnail 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: 104 (solid), 105 (dashed), 106 (dotted), and 107 (dot-dashed).

For the models with a Salpeter mass function, we have also evaluated the mass-dependent 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 104tdyn, 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 multi-particle 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 Osipkov-Merritt 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 rs, or of about five final half-mass radii r50%.

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 (MIMBH/⟨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−1r50%. A result whose importance cannot be overstated is that the distribution of the distance as well as the sky-projected 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.


1

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).

2

We note that, in Osipkov-Merritt models, the radial profile of β can be written explicitly as a function of the anisotropy radius ra as .

3

Such a radius of about 0.1rs is also consistent with the typical core-stalling 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 MIUR-PRIN2017 project Coarse-grained description for non-equilibrium systems and transport phenomena (CO-NEST) 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 Sklodowska-Curie grant agreement No. 896248. We warmly acknowledge Evangelia Tremou for reading a draft of this manuscript.

References

  1. Aros, F. I., Sippel, A. C., Mastrobuono-Battisti, A., et al. 2020, MNRAS, 499, 4646 [CrossRef] [Google Scholar]
  2. Bahcall, J. N., & Wolf, R. A. 1976, ApJ, 209, 214 [NASA ADS] [CrossRef] [Google Scholar]
  3. Ballone, A., Mapelli, M., & Pasquato, M. 2018, MNRAS, 480, 4684 [CrossRef] [Google Scholar]
  4. Banik, U., & van den Bosch, F. C. 2021, ApJ, 912, 43 [NASA ADS] [CrossRef] [Google Scholar]
  5. Banik, U., & van den Bosch, F. C. 2022, ApJ, 926, 215 [NASA ADS] [CrossRef] [Google Scholar]
  6. Baumgardt, H. 2001, MNRAS, 325, 1323 [NASA ADS] [CrossRef] [Google Scholar]
  7. Baumgardt, H., & Hilker, M. 2018, MNRAS, 478, 1520 [Google Scholar]
  8. Baumgardt, H., De Marchi, G., & Kroupa, P. 2008, ApJ, 685, 247 [NASA ADS] [CrossRef] [Google Scholar]
  9. Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition (Princeton University Press) [Google Scholar]
  10. Bortolas, E., Gualandris, A., Dotti, M., Spera, M., & Mapelli, M. 2016, MNRAS, 461, 1023 [Google Scholar]
  11. Breen, P. G., Varri, A. L., & Heggie, D. C. 2017, MNRAS, 471, 2778 [NASA ADS] [CrossRef] [Google Scholar]
  12. Brockamp, M., Baumgardt, H., & Kroupa, P. 2011, MNRAS, 418, 1308 [CrossRef] [Google Scholar]
  13. Bustamante-Rosell, M. J., Noyola, E., Gebhardt, K., et al. 2021, ApJ, 921, 107 [NASA ADS] [CrossRef] [Google Scholar]
  14. Carr, B. J., Bond, J. R., & Arnett, W. D. 1984, ApJ, 277, 445 [NASA ADS] [CrossRef] [Google Scholar]
  15. Chandrasekhar, S., & von Neumann, J. 1942, ApJ, 95, 489 [NASA ADS] [CrossRef] [Google Scholar]
  16. Chandrasekhar, S., & von Neumann, J. 1943, ApJ, 97, 1 [NASA ADS] [CrossRef] [Google Scholar]
  17. Chatterjee, P., Hernquist, L., & Loeb, A. 2002a, Phys. Rev. Lett., 88, 121103 [NASA ADS] [CrossRef] [Google Scholar]
  18. Chatterjee, P., Hernquist, L., & Loeb, A. 2002b, ApJ, 572, 371 [NASA ADS] [CrossRef] [Google Scholar]
  19. Ciotti, L. 1996, ApJ, 471, 68 [NASA ADS] [CrossRef] [Google Scholar]
  20. Dehnen, W., & Read, J. I. 2011, Eur. Phys. J. Plus, 126, 55 [Google Scholar]
  21. Dejonghe, H. 1987, MNRAS, 224, 13 [NASA ADS] [CrossRef] [Google Scholar]
  22. Di Carlo, U. N., Mapelli, M., Pasquato, M., et al. 2021, MNRAS, 507, 5132 [NASA ADS] [CrossRef] [Google Scholar]
  23. Di Cintio, P., Livi, R., Lepri, S., & Ciraolo, G. 2017, Phys. Rev. E, 95, 043203 [NASA ADS] [CrossRef] [Google Scholar]
  24. 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]
  25. Di Cintio, P., Pasquato, M., Kim, H., & Yoon, S.-J. 2021, A&A, 649, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Di Cintio, P., Pasquato, M., Simon-Petit, A., & Yoon, S.-J. 2022, A&A, 659, A19 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Di Cintio, P., Pasquato, M., Barbieri, L., et al. 2023, IAU Symp., 362, 134 [NASA ADS] [Google Scholar]
  28. Ebisuzaki, T., Makino, J., Tsuru, T. G., et al. 2001, ApJ, 562, L19 [NASA ADS] [CrossRef] [Google Scholar]
  29. Freitag, M., & Benz, W. 2001, A&A, 375, 711 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  30. Gieles, M., & Zocchi, A. 2015, MNRAS, 454, 576 [NASA ADS] [CrossRef] [Google Scholar]
  31. Giersz, M. 2006, MNRAS, 371, 484 [NASA ADS] [CrossRef] [Google Scholar]
  32. Giersz, M., Heggie, D. C., Hurley, J. R., & Hypki, A. 2013, MNRAS, 431, 2184 [NASA ADS] [CrossRef] [Google Scholar]
  33. Greene, J. E., Strader, J., & Ho, L. C. 2020, ARA&A, 58, 257 [Google Scholar]
  34. Greif, T. H. 2015, Comput. Astrophys. Cosmol., 2, 3 [NASA ADS] [CrossRef] [Google Scholar]
  35. Heggie, D. C. 2011, Bull. Astron. Soc. India, 39, 69 [NASA ADS] [Google Scholar]
  36. Holley-Bockelmann, K., Gültekin, K., Shoemaker, D., & Yunes, N. 2008, ApJ, 686, 829 [Google Scholar]
  37. Holtsmark, J. 1919, Ann. Phys., 363, 577 [CrossRef] [Google Scholar]
  38. Hurley, J. R., Pols, O. R., Aarseth, S. J., & Tout, C. A. 2005, MNRAS, 363, 293 [NASA ADS] [CrossRef] [Google Scholar]
  39. Hypki, A., & Giersz, M. 2013, MNRAS, 429, 1221 [NASA ADS] [CrossRef] [Google Scholar]
  40. Inayoshi, K., Visbal, E., & Haiman, Z. 2020, ARA&A, 58, 27 [NASA ADS] [CrossRef] [Google Scholar]
  41. Kamlah, A. W. H., Leveque, A., Spurzem, R., et al. 2022, MNRAS, 511, 4060 [NASA ADS] [CrossRef] [Google Scholar]
  42. Kim, E., Yoon, I., Lee, H. M., & Spurzem, R. 2008, MNRAS, 383, 2 [NASA ADS] [CrossRef] [Google Scholar]
  43. Li, G.-P. 2022, A&A, 666, A194 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  44. Lodato, G., & Natarajan, P. 2006, MNRAS, 371, 1813 [NASA ADS] [CrossRef] [Google Scholar]
  45. Loeb, A., & Rasio, F. A. 1994, ApJ, 432, 52 [Google Scholar]
  46. Londrillo, P., & Messina, A. 1990, MNRAS, 242, 595 [NASA ADS] [CrossRef] [Google Scholar]
  47. Merritt, D. 1985, AJ, 90, 1027 [Google Scholar]
  48. Merritt, D. 2004, in Coevolution of Black Holes and Galaxies, ed. L. C. Ho, 263 [Google Scholar]
  49. Mezcua, M. 2017, Int. J. Mod. Phys. D, 26, 1730021 [Google Scholar]
  50. Miller, M. C., & Hamilton, D. P. 2002, MNRAS, 330, 232 [Google Scholar]
  51. Osipkov, L. P. 1979, Soviet. Astron. Lett., 5, 42 [NASA ADS] [Google Scholar]
  52. Pacucci, F., & Loeb, A. 2022, MNRAS, 509, 1885 [Google Scholar]
  53. Pasquato, M., & Di Cintio, P. 2020, A&A, 640, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  54. Peebles, P. J. E. 1972, ApJ, 178, 371 [NASA ADS] [CrossRef] [Google Scholar]
  55. Plummer, H. C. 1911, MNRAS, 71, 460 [Google Scholar]
  56. Portegies Zwart, S. F., Baumgardt, H., Hut, P., Makino, J., & McMillan, S. L. W. 2004, Nature, 428, 724 [Google Scholar]
  57. Rizzuto, F. P., Naab, T., Spurzem, R., et al. 2021, MNRAS, 501, 5257 [NASA ADS] [CrossRef] [Google Scholar]
  58. Rizzuto, F. P., Naab, T., Rantala, A., et al. 2023, MNRAS, 521, 2930 [NASA ADS] [CrossRef] [Google Scholar]
  59. Salpeter, E. E. 1955, ApJ, 121, 161 [Google Scholar]
  60. Sollima, A., & Ferraro, F. R. 2019, MNRAS, 483, 1523 [Google Scholar]
  61. Sollima, A., & Mastrobuono Battisti, A. 2014, MNRAS, 443, 3513 [NASA ADS] [CrossRef] [Google Scholar]
  62. Takahashi, K., & Portegies Zwart, S. F. 2000, ApJ, 535, 759 [NASA ADS] [CrossRef] [Google Scholar]
  63. Torniamenti, S., Pasquato, M., Di Cintio, P., et al. 2022, MNRAS, 510, 2097 [Google Scholar]
  64. Tremou, E., Strader, J., Chomiuk, L., et al. 2018, ApJ, 862, 16 [Google Scholar]
  65. Vasiliev, E. 2015, MNRAS, 446, 3150 [NASA ADS] [CrossRef] [Google Scholar]
  66. Volonteri, M., Lodato, G., & Natarajan, P. 2008, MNRAS, 383, 1079 [Google Scholar]
  67. Wang, L. 2020, MNRAS, 491, 2413 [NASA ADS] [Google Scholar]
  68. Weatherford, N. C., Kıroğlu, F., Fragione, G., et al. 2023, ApJ, 946, 104 [NASA ADS] [CrossRef] [Google Scholar]
  69. Webb, J. J., Hunt, J. A. S., & Bovy, J. 2023, MNRAS, 521, 3898 [NASA ADS] [CrossRef] [Google Scholar]
  70. Yoshida, N., Omukai, K., Hernquist, L., & Abel, T. 2006, ApJ, 652, 6 [NASA ADS] [CrossRef] [Google Scholar]
  71. Zocchi, A., Gieles, M., & Hénault-Brunet, 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

Table 1.

Summary of the initial conditions.

Table 2.

Summary of the simulation properties.

All Figures

thumbnail Fig. 1.

Observational claims of IMBH detection in the (μ = MIMBH/⟨m⟩, α = MIMBH/M) plane, shown as black diamonds. Any technique (most notably the direct N-body 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 N-body codes currently available, N = 106. As a reference, N = 105 is also shown as a dashed line.

In the text
thumbnail Fig. 2.

Evolution of the Lagrangian radii r5% enclosing 5% of the system’s mass (top panels) and central velocity dispersion σc evaluated for particles inside r5% (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 = 105 so that μ = α × 105.

In the text
thumbnail Fig. 3.

Evolution of the radial anisotropy profile β(r) in Plummer models without a central IMBH, with a Salpeter IMF, N = 106, 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 dot-dashed lines mark the profile β at the indicated time of core collapse tcc.

In the text
thumbnail Fig. 4.

Evolution of the radial anisotropy profile β(r) in Plummer models with a central IMBH of mass MIMBH = 10−3M, with a Salpeter IMF (top panel row) or equal masses (bottom panel row), N = 106, 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
thumbnail 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 = 105 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
thumbnail 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 = 105.

In the text
thumbnail Fig. 7.

Distributions of the radial coordinate (left panel) and velocity (middle panel) of the IMBHs embedded in a Plummer model with 3 × 103 ≤ N ≤ 106 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 rwan has a trend compatible with μ−1/2 for μ < 3 × 102 (dash ed line). In all cases, α = 10−3.

In the text
thumbnail Fig. 8.

Wander radius (scaled to 0-1 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: 104 (solid), 105 (dashed), 106 (dotted), and 107 (dot-dashed).

In the text
thumbnail Fig. 9.

Distribution of the IMBH radial coordinate for Plummer models with N = 106, MIMBH = 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
thumbnail 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: 104 (solid), 105 (dashed), 106 (dotted), and 107 (dot-dashed).

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.