Intermediate-Mass Black Holes in binary-rich star clusters

TL;DR

This study uses N-body simulations to explore how intermediate-mass black holes influence the dynamics of binary-rich globular clusters, revealing signatures like high-velocity escapers and increased star-IMBH interactions.

Contribution

It demonstrates that IMBHs create observable dynamical effects in binary-rich clusters, including high-velocity stars and enhanced star-IMBH interactions, advancing methods to detect IMBHs.

Findings

IMBHs develop a star density cusp around them.

Presence of binaries leads to high-velocity star escapers.

Binary populations increase star-IMBH interaction rates.

Abstract

There is both theoretical expectation and some observational clues that intermediate mass black holes reside in nuclei of globular clusters. In order to find an independent indicator for their existence, we investigate in this paper how an IMBH manifests itself through its dynamical interaction with a binary rich globular cluster of moderate extension and mass. By means of direct $N$-body integration we follow the dynamical evolution of models of such a system over a time span of $\approx$ 0.8 Gyr and compare the cases with and without the primordial binaries as well as with and without the IMBH. In accord with previous results, we show that when present the IMBH develops a power-law density cusp of stars around it, regardless of the binary population in the cluster. If, however, binaries are present, their interaction with the IMBH leads to the production of high velocity escapers at a rate of the order of 0.1 Myr$^{-1}$. These stars may contribute to the population of high-velocity stars observed in the Galaxy. Clusters hosting the IMBH together with high number of binaries also form a denser halo of marginally unbound stars than clusters that lack either the IMBH or the rich binary population. Finally, we show that the binary population leads to an increased rate of direct interactions of stars with the IMBH, potentially observable as tidal disruption events.

Tables

Table 1. Table 1: Variable parameters of the models: name, number of stars ( N 𝑁 N ), mass of the IMBH ( M IMBH subscript 𝑀 IMBH M_{\mathrm{IMBH}} ) and primordial binary fraction ( f bin subscript 𝑓 bin f_{\mathrm{bin}} ).

name	$N$	$M_{\mathrm{IMBH}}$	$f_{\mathrm{bin}}$
BH-BIN	$50001$	${10}^3{M}_{\odot }$	$0.5$
BH	$50001$	${10}^3{M}_{\odot }$	$0.0$
BIN	$50000$	$0$	$0.5$

Equations

$$r_i\propto \frac{M_{\mathrm{IMBH}}}{M_{\mathrm{GC}}}{r}_c.$$

$$r_t\approx {\left(\frac{M_{\mathrm{IMBH}}}{m}\right)}^{1/3}a.$$

$$v_{ej}\approx 460{\left(\frac{a}{0.1\mathrm{AU}}\right)}^{-1/2}{\left(\frac{m}{2{\mathrm{M}}_{\odot }}\right)}^{1/3}{\left(\frac{M_{\mathrm{IMBH}}}{10^3{\mathrm{M}}_{\odot }}\right)}^{1/6}\mathrm{km}{\mathrm{s}}^{-1}.$$

$$v_1=v_{ej}{\left(\frac{2m_2}{m}\right)}^{1/2}\text{makebox}[30.00005pt]and{v}_2=v_{ej}{\left(\frac{2m_1}{m}\right)}^{1/2},$$

$$\mathcal{R}$$

$$\mathcal{R}\approx f_b{\left(\frac{n}{10^5{\mathrm{pc}}^3}\right)}^2{\left(\frac{M_{\mathrm{IMBH}}}{10^3{\mathrm{M}}_{\odot }}\right)}^3{\left(\frac{\sigma }{10\mathrm{km}{\mathrm{s}}^{-1}}\right)}^{-9}{\mathrm{Myr}}^{-1}.$$

$$\frac{\mathrm{d}{N}_{\mathrm{esc}}}{\mathrm{d}{v}_{\mathrm{esc}}\mathrm{d}t}\approx 350{\left(\frac{v_{\mathrm{esc}}}{\mathrm{km}{\mathrm{s}}^{-1}}\right)}^{-3}{\left(\mathrm{km}{\mathrm{s}}^{-1}\right)}^{-1}{\mathrm{Myr}}^{-1}.$$

$${\Phi }_{\mathrm{BH}}(r)=-\frac{G{M}_{\mathrm{BH}}}{r},$$

$${\Phi }_{\mathrm{b}}(r)=-\frac{G{M}_{\mathrm{bul}}}{r+a},$$

$${\Phi }_{\mathrm{disk}}(R,z)=-\frac{G{M}_{\mathrm{disk}}}{\sqrt{(}{R}^2+(b+\sqrt{c^2+z^2}{)}^2)},$$

$${\Phi }_{\mathrm{halo}}(r)=-\frac{G{M}_{\mathrm{DM}}\ln (1+r/r_s)}{r}.$$

$$2v_{\mathrm{esc}}^2(r)=-V(r)=\frac{G}{r}\sum _{j=1}^i{m}_j+\sum _{j=i+1}^N\frac{G{m}_j}{r_j},$$

Full text

1 Introduction

Independently of their spin (and possibly charge), black holes are usually divided into three different categories. Supermassive black holes (SMBHs, $M\gtrsim 10^{5}\,\rm M_{\odot}$) reside in the centres of galaxies, where they shape the surrounding gas and stellar content (Hopman & Alexander, 2006; Alexander, 2017). Stellar-mass black holes (SBHs, $10\,\rm M_{\odot}\lesssim M\lesssim 100\,\rm M_{\odot}$) are the final stage of massive stars lifetime and have been recently observed by LIGO detections in SBH-SBH merger events (Abbott et al., 2016; Abbott et al., 2017; Fragione & Kocsis, 2018). However, for the intermediate-size population, the so-called intermediate-mass black holes (IMBHs) ($100\,\rm M_{\odot}\lesssim M\lesssim 10^{5}\,\rm M_{\odot}$), there is not yet a definitive accepted proof of their existence. In the Milky Way, two clusters have been claimed to host an IMBH based on dynamical modeling, namely $\omega$ Cen (Noyola et al., 2010; Jalali et al., 2012; Baumgardt, 2017) and 47 Tuc (Kızıltan et al., 2017), but other authors (Anderson & van der Marel, 2010; van der Marel & Anderson, 2010; Zocchi et al., 2017, 2018) have advised caution with this interpretation. However, the recent observation of a tidal disruption event (TDE) consistent with an IMBH of $\sim 5\times 10^{4}\,\rm M_{\odot}$ in an off-centre star cluster, at a distance of $\sim 12.5$ kpc from the centre of the host galaxy, has provided a new, maybe final, proof of their existence (Lin et al., 2018).

Assuming that the observed $M_{\mathrm{SMBH}}-\sigma$ relation ($\sigma$ is the velocity dispersion) holds also in the range of IMBH masses (Merritt, 2013), the ideal place to look for IMBHs is at the centres of globular clusters (GCs). In GCs, the most massive stars may segregate and merge in the core, forming an object in the mass range of an IMBH (Freitag et al., 2006). Portegies Zwart & McMillan (2002) suggested that star clusters with small initial half-mass relaxation times may form an IMBH with a mass up to $\sim 0.1$% of the mass of the entire star cluster, as a consequence of the high stellar collision rates. Recently, Giersz et al. (2015) found slow and fast scenarios of IMBH mass growth: GCs with a large initial cluster concentration have a larger probability of forming an IMBH earlier and faster. Once formed, an IMBH resides at the centre of its host GC (Chatterjee et al., 2002a, b), and interacts with the host cluster stars and compact objects (Leigh et al., 2014). Unfortunately, both the lack of high-resolution data and accurate N-body modeling of star clusters deny definitive kinematic detections of IMBH in the centres of GCs.

Galactic nuclei may host IMBHs as well. The centre of the Milky Way may host an IMBH in the central parsec and several IMBHs in its nuclear star cluster, possibly delivered by inspiraling star clusters (Portegies Zwart & McMillan, 2002; Fragione, Ginsburg & Kocsis, 2018; Fragione, Leigh, Ginsburg & Kocsis, 2018). In galactic nuclei in general, the TDE rate may be enhanced due to the presence of an IMBH (Chen et al., 2009, 2011; Fragione & Leigh, 2018b).

IMBHs may be illuminated by Gravitational Wave (GW) events. IMBH-SBH binaries may commonly form in the core of GCs and may merge as intermediate mass ratio inspirals (IMRIs), at a rate as large as $\sim$ few $\mathrm{Gpc}^{-3}\,\mathrm{yr}^{-1}$. IMBH-SBH binaries of different masses can be observed by present and upcoming facilities, such as LIGO, the Einstein Telescope, and LISA (Leigh, Lützgendorf, Geller, Maccarone, Heinke & Sesana, 2014; Amaro-Seoane & Santamaría, 2010; Fragione, Ginsburg & Kocsis, 2018; Fragione & Leigh, 2018a). Besides GW events, TDEs may also take place in the center of a cluster hosting an IMBH (Baumgardt, Makino & Ebisuzaki, 2004a, b; Fragione, Leigh, Ginsburg & Kocsis, 2018), whose rate ($10^{-6}$-$10^{-3}$ yr*-1*) can be of the same order of the TDEs driven by the host galaxy SMBH.

Several efforts have been made in modeling GCs with IMBHs through direct $N$-body simulations. Baumgardt et al. (2005) investigated the radial density profile of stellar systems with a central IMBH. Trenti et al. (2007) and Gill et al. (2008) used direct N-body simulations to study clusters with an IMBH and $5$-$10$% primordial binaries. One of the predictions of these models was that clusters hosting an IMBH can develop a cusp of stars in the projected surface brightness profile (Baumgardt et al., 2005; Miocchi, 2007), which is usually regarded as an indicator of the presence of an IMBH, even though also other dynamical processes can produce similar signatures (see e.g. Hurley, 2007; Vesperini et al., 2010). Lützgendorf et al. (2013) performed $N$-body simulations of GCs hosting IMBHs as a function of the neutron star and SBH retention fraction, while Leigh et al. (2014) studied the co-existence of SBH binaries with an IMBH. Recently, Baumgardt (2017) has run a large grid of N-body simulations with different mass fractions $M_{\mathrm{IMBH}}/M_{\mathrm{GC}}$ (where $M_{\mathrm{IMBH}}$ and $M_{\mathrm{GC}}$ are the IMBH and GC mass, respectively).

Presence of large number of primordial binaries affects evolution of star clusters, namely their cores during the collapse and subsequent gravo-thermal oscillations (see, e.g., Heggie et al., 2006; Trenti et al., 2007, for extensive study on this topic). It is also expected that interaction of binaries with an IMBH leads to additional phenomena, e.g., ejection of high velocity stars through the Hills mechanism, similarly to what happens in Galactic Centre (Hills, 1988; Yu & Tremaine, 2003; Bromley et al., 2006; Sari et al., 2010; Brown et al., 2018). In this scenario, a binary star passing close to the IMBH is tidally separated with one component being ejected away at high velocity, while the other one remains tightly bound to the IMBH (Fragione & Gualandris, 2018a).

A full $N$-body treatment of star cluster with primordial binaries and the IMBH is a computationally challenging task and, therefore, only few works on this topic can be found in the literature so far (see, e.g., Trenti et al., 2007, for one of the first attempts). In this paper, we aim to further fill in this gap, presenting the first direct $N$-body simulations of star clusters hosting an IMBH and a large primordial binary fraction (up to a binary fraction of $50$%). We use for our simulations a modified version of the NBODY6 code (Aarseth, 2003). We study clusters with or without an IMBH and with different binary fractions, show how the presence of the IMBH affects the central parts of the cluster and discuss the rate at which the IMBH ejects stars out of the cluster into the host galaxy. We then discuss how these ejected stars may contribute to the observed population of hypervelocity (HVSs) and runaway stars (RSs), along with the other mechanisms proposed in literature (Hills, 1988; Yu & Tremaine, 2003; Gualandris & Portegies Zwart, 2007; Sesana et al., 2006; Zubovas et al., 2013; Capuzzo-Dolcetta & Fragione, 2015; Fragione & Capuzzo-Dolcetta, 2016; Šubr & Haas, 2016; Boubert et al., 2017; Fragione et al., 2017; Fragione & Gualandris, 2018b, a).

The paper is organized as follows. In Section 2, we present the numerical method we use to evolve the GC. In Section 3, we present our results. Finally, in Section 4, we draw our conclusions.

2 Method

Our numerical simulations with the NBODY6 code (Aarseth, 2003) aim at studying the impact of an IMBH on the properties and the evolution of an aged stellar system like a low-mass GC. The initial structure, stellar populations and binary populations are therefore chosen to represent roughly such a system, while avoiding too many speculative details and keeping the models numerically manageable. IMBHs are thought to form early in the evolution of a massive star cluster (see Section 1), if they form at all, and their formation is therefore not under scrutiny in the present models. The models are rather set up with a particle representing the IMBH, and compared to models that are set up without an IMBH, but identical parameters otherwise. This allows to identify the impact of the IMBH on such systems, and how the presence of an IMBH may be detected observationally.

The simultaneous treatment of a large population of binaries, as it is likely to be found in GCs (like in any known stellar population), and an IMBH has not been done in previous simulations. The reason is that binary systems are numerically demanding and large differences in the masses of the simulated particles (as they are implied by the presence of an IMBH among stars) are a challenge for direct $N$-body integrators. The simultaneous treatment of these two complications is what limits our simulations to 50000 stars, given that an additional constraint is that for meaningful results on the long-term evolution and stability of these systems, we need to simulate several hundreds of Myr of cluster evolution within reasonable computing time. Thus, 50000 particles sounds like a modest number in comparison to the $10^{6}$ particles considered by Wang et al. (2016), but we note that in the simulations of Wang et al. (2016) only $5$% of the stars are in primordial binaries (in contrast to $50$% in our simulations), and also the IMBH is missing.

2.1 Numerical model for a star cluster with binaries and an IMBH

In the following, we give a detailed description of our setup for the system that is at the centre of attention in the present paper, namely a massive intermediate-age star cluster with a large population of binaries and an IMBH. We will refer to this model as model BH-BIN, and name and discuss our choices regarding the stellar mass function (Section 2.1.1), the initial binary population (Section 2.1.2), and the initial cluster structure (Section 2.1.3) in this model. In Section 2.2, we will briefly describe some additional $N$-body models that we run in order to evaluate which impact the inclusion of a large binary population and an IMBH has on the evolution of a star cluster.

2.1.1 Mass function

We assume for model BH-BIN that the initial masses of stars are distributed according to a broken power-law distribution function with the power-law index $\alpha=-1.3$ for $M_{\star}<0.5\,\rm M_{\odot}$ and $\alpha=-2.3$ for $0.5\,\rm M_{\odot}\leq M_{\star}<m_{\rm max}$. This mass function is known as the canonical stellar initial mass function, and is characteristic for many star formation events, except the most extreme ones (Kroupa, 2001; Kroupa et al., 2013).

Furthermore, all stars are assumed to have formed instantaneously in a single event. In reality, stars certainly do not form instantaneously in a star cluster, and some models for explaining the observed element abundances in stars in GCs even imply that the stars in GCs were born in multiple star formation epochs instead of a single one. However, these events take place on a time scale of several 100 Myrs at most, so that the birth of all stars in an initial star burst becomes a very good approximation when the stellar masses of a several Gyr old GC are considered.

The initial number of particles representing individual stars and stellar remnants is 50000. Of these particles, 44500 represent stars. Their masses are drawn randomly from the canonical IMF in the mass interval $(0.1\ {\rm M}_{\odot},\,1\ {\rm M}_{\odot})$. The upper limit of assumed stellar masses is motivated with the finding that already in a 1 Gyr old stellar population, the most massive surviving stars have a mass below $1\,\,\rm M_{\odot}$, and from then onwards, the upper mass limit for surviving stars keeps on dropping only slowly as the stellar population continues to age (see e.g. the simple stellar population models by Maraston 2005). Our upper mass limit therefore roughly represents the situation in a few Gyr old star cluster.

The remaining 5500 particles with stellar masses represent stellar remnants; i.e. white dwarfs, neutron stars and stellar-mass black holes. Their number corresponds to the number of stars with a mass above $1\,\,\rm M_{\odot}$ according to the canonical IMF, if the number of stars below $1\,\,\rm M_{\odot}$ is 44500, and the upper mass limit of the IMF, $m_{\rm max}$ is set to $100\,\,\rm M_{\odot}$. The actual upper mass limit of the IMF may be significantly higher (cf. Crowther et al. 2010), but if the IMF is canonical, stars with initial masses above 100 $\,\rm M_{\odot}$ are very rare and can therefore be neglected for the purpose of this paper.

The mass attributed to each particle representing a stellar remnant is $1\,\,\rm M_{\odot}$. This is only an approximation, but it seems justifiable, given that the observed masses of white dwarfs in Kalirai et al. (2008) range between $0.5\,\,\rm M_{\odot}$ and $1.1\,\,\rm M_{\odot}$, and given that the observed masses of neutron stars seem all close to a value of $1.35\,\,\rm M_{\odot}$ (Thorsett & Chakrabarty, 1999). Stars with an initial mass above $25\,\,\rm M_{\odot}$ are likely to become stellar-mass black holes which can be much heavier than $1\,\,\rm M_{\odot}$ (compare for example figures 12 and 16 in Woosley et al. 2002), but their numbers are suppressed strongly by the steep slope of the canonical IMF at high masses.

We finally include also one additional particle with a mass of $1000\,\,\rm M_{\odot}$, which represents the IMBH. This additional particle is initially placed at a random position in the cluster, but sinks rapidly towards the centre of the cluster through dynamical friction due to its high mass (Binney & Tremaine, 2008).

Thus, in essence, the model BH-BIN is set up to represent few Gyr old massive star clusters in which stars were formed according to the canonical IMF. Stellar evolution is considered by the choice of the upper mass limit for still luminous stars and the ratio between the number of luminous stars and stellar remnants, while dynamical evolution, which leads to a loss of low-mass stars (Baumgardt & Makino, 2003), is neglected.

2.1.2 Binaries

Kroupa (1995) argues that initially nearly every star in a star cluster is part of a binary or a multiple system of even higher order. However, many of these initial binaries are very wide and therefore easily disrupted through encounters. Thus, in dynamically evolved star clusters, the fraction of stars in multiple systems will be much lower.

The actual properties of the binary populations in GCs in not known in much detail, but as an approximation to the binary population in GCs, another population that has undergone strong dynamical processing may serve. If essentially all stars were born in star clusters, as suggested by Lada & Lada (2003), such a population is found in the Galactic field, because the stars in the Galactic field then all come from star clusters that were dissolved through dynamical processes. The properties of binaries in the Galactic field are well studied by Duquennoy & Mayor (1991), which is therefore chosen as the basis for the present work. We note that according to Minor (2013) also the properties of the binary populations in the dwarf Spheroidal Galaxies which he studied are consistent with those Duquennoy & Mayor (1991) found for the Galactic field.

The binary fraction, $f_{\mathrm{bin}}$, is defined as the number of binaries over the number of center-of-mass systems, i.e. binaries and single stars in our case. We set $f_{\mathrm{bin}}=0.5$ in our model BH-BIN, which is close to the value Duquennoy & Mayor (1991) find in the Galactic field (but much lower than $f_{\mathrm{bin}}=1$, which was suggested by Kroupa (1995) for dynamically unevolved star clusters). We note that Wang et al. (2016) assumed $f_{\mathrm{bin}}=0.05$ in their N-Body simulation with $10^{6}$ individual stars, while Sollima et al. (2007) found $f_{\mathrm{bin}}>0.06$ for the low-density GCs they studied, and in some cases $f_{\mathrm{bin}}\approx 0.5$ seems likely. This motivates our high choice for $f_{\mathrm{bin}}$ in model BH-BIN. However, in more massive and thus (for a fixed radius) more dense GCs, the binary fraction appears to be smaller ($0.04<f_{\mathrm{bin}}<0.09$, see Cool & Bolton 2002 for the case of NGC 6397 and Romani & Weinberg 1991 for the case of M92; see also Sollima et al. 2007), and thus the choice of Wang et al. (2016) seems appropriate for the model they consider.

The semi-major axes of primordial binaries follow the distribution found by Duquennoy & Mayor (1991) for the Galactic field in all models we consider, i.e. a log-normal distribution which peaks at a semi-major axis corresponding to an orbital period of $\sim 180$ yr. The interpretation of this distribution according to Kroupa (1995) is that binaries with very wide orbits, which may have been frequent initially, are already destroyed through dynamical interactions. The initial eccentricities of the binaries are set to zero for simplicity.

The procedure of pairing is based on random number generator. All objects (stars and compact remnants) are first put into a pool ordered according to their mass from the most massive one to the lightest. Starting from the most massive object, we generate a random number which determines whether it will be paired with another one or not (with equal probability for these two cases). If the object is to be a member of a binary, its secondary is selected with uniform probability from the objects of the same type remaining in the pool (i.e., binaries made of main sequence star and compact remnant are not allowed). Both objects are then removed from the pool and the procedure continues on, until all stars are taken from the pool. We do not set up higher-order systems.

2.1.3 Initial cluster structure

In our setup for model BH-BIN, the star cluster initially follows the Plummer density profile with half-mass radius, $r_{\mathrm{h}}=3{\,\rm pc}$. The Plummer density profile has been shown by Plummer (1911) to provide simple, yet very satisfactory fits to the observed density profiles of GCs, and is therefore still a very popular choice for setting up models for star clusters.

The chosen half-mass radius is quite characteristic for GCs in the Milky Way (see figure 8 in McLaughlin 2000), but also for extragalactic GCs in the Virgo Cluster (Jordán et al., 2005) and the Fornax Cluster (Jordán et al., 2015). The choice for the half-mass radius puts our model, together with the mass attributed by 50000 stars (see Section 2.1.1), in the regime of low-mass and low-density GCs. This makes the choice of rather high binary fractions suggested in Sollima et al. (2007) for low-density appropriate for them, while the case of a low binary fraction that Wang et al. (2016) use in their models of rather massive GCs is depreciated for our case (see also Section 2.1.2).

The initial positions and velocities were generated using the plumix code (Šubr et al., 2008; Šubr, 2012) which ensures that the IMBH local mass overdensity is somewhat reduced by lower local density of lighter objects.

Finally, all the cluster models are generated in isolation, i.e. without any external tidal field.

2.2 Other N-body models

The focus of this work lies on the evolution of massive star clusters with an IMBH and binaries, which is the situation depicted with our model BH-BIN. However, in order to evaluate how the interplay of an IMBH and a large population of primordial binaries affects the evolution of a massive star cluster, we have integrated additional N-body models, which we list below:

• •

Model BIN is the same as model BH-BIN, except that there is no particle representing the IMBH included in model BIN.

• •

Model BH is the same as model BH-BIN, except that $f_{\rm bin}=0$ in model BH, so that there are initially no binaries in the simulation, while they may form dynamically later on.

A summary and a quick reference to our models is given in Table (1). For each model, seven realisations which differ in initial positions and velocities of stars, but sharing the overall statistical properties, were integrated. Results presented below are based on averages over the seven realisations of individual models.

2.3 Numerical integrator

We used NBODY6 code (Aarseth, 2003) for the numerical integration of the equations of motion. We added to the original code routines for logging of beginnings and endings of regularizations into a binary file. We also altered the decision making algorithm for adding particles to neighbor lists. In particular, we weight the standard distance criterion by the mass of the given particle so that the more massive particles are added to the list even when they are at larger distances than the lighter ones. The most prominent target was the IMBH particle which was, due to this modification, a member of the neighbor lists of all other (star) particles. The modification of the neighbor list influences the integration and increases its stability. Among many runtime options of the NBODY6 code, let us specifically mention that we switched off the internal evolution of the stars as well as the post-Newtonian corrections to the stellar dynamics.

3 Results

3.1 Radial density profile

The IMBH which was set at random position in the star cluster initially, sinks quickly into the cluster centre on the time scale $\lesssim 20{\,\rm Myr}$ due to dynamical friction (Binney & Tremaine, 2008). From that time on, a stellar cusp evolves around the IMBH and it is clearly visible in the radial density profile already at $T\geq 100{\,\rm Myr}$ (see Fig. 1). The buildup of the cusp is somewhat faster in model BH-BIN than in BH, nevertheless, already at $T\gtrsim 200{\,\rm Myr}$, the difference between the two models becomes negligible. We also note that the region below $\approx 0.1{\,\rm pc}$ is occupied by no more than several tens of stars and, therefore, the radial density there fluctuates in time and individual realisations of a particular model may differ from each other. Once the central cusp approaches the equilibriums state, its evolution slows down and both models that include the IMBH share a similar density profile. The central cusp (within the IMBH sphere of influence $r_{i}$) is well approximated by the Bahcall & Wolf (1976) density profile $\varrho(r)\propto r^{-7/4}$. In agreement with previous works (Baumgardt et al., 2005; Miocchi, 2007), we find that a shallow cusp in the projected surface brightness profile is a possible indicator of the presence of an IMBH, but other dynamical processes can produce cuspy profiles as well (Hurley, 2007; Vesperini et al., 2010).

The fact that the radial density profile reaches near equilibrium state in the innermost parts of the models including the IMBH, however, does not mean that the central cusp does not lead to any interesting phenomena. At each moment after $T\approx 100{\,\rm Myr}$, there are typically one or two stars or stellar remnants orbiting the IMBH with semi-major axis smaller than $\lesssim 10^{-4}{\,\rm pc}$ and eccentricity $e>0.99$ (see also MacLeod et al., 2016, for discussion of dynamical formation and properties of close stellar companions to an IMBH). Peak velocities of these stars during passages through the percentres of their orbits are few thousands of kilometers per second, i.e., comparable to velocities of so-called S-stars orbiting the Galactic SMBH (see, e.g., Ghez et al., 2005; Gillessen et al., 2009). Integration of these stars is in general time demanding due to large accelerations and the fact that the tightly bound system consisting of more than two particles (including the IMBH) is not suitable for standard regularisation techniques implemented in $N$-body integrators. Despite their large binding energy, even the orbits of the most tightly bound stars are being subject to perturbations which often lead to their collision with the IMBH.

The region at distances $10^{-4}{\,\rm pc}\lesssim r\lesssim 0.1{\,\rm pc}$ from the IMBH is typically occupied by several tens of solar-mass objects which can be considered bound to the IMBH. Nevertheless, it appears that the time an individual star spends on the orbit bound to the IMBH is quite short – more than 60% of stars are scattered out from this region on the time-scale of $\approx 16{\,\rm Myr}$. This quite strong relaxation in the vicinity of the IMBH is very likely another source of computational error which leads to slow down of the integrations.

In our models, we have adopted an initial Plummer profile for the initial density profile of the clusters. Since the cluster loses memory of the initial conditions after roughly one relaxation time, the long-term evolution of the cluster does not depend significantly on the details of the initial conditions, once the IMBH mass and the primordial binary fractions are fixed. Trenti et al. (2007) found that concentrated King models initially presents a core expansion due to the high energy injected by three- and four-body interactions, while clusters with a shallow King profile or a Plummer profile have an initial contraction. Nevertheless, both concentrated and shallow profiles tend to a common value of the core radius $r_{c}$. Typically, the size of the IMBH sphere of influence is proportional to the core radius of the host GC (Baumgardt et al., 2004a, b)

[TABLE]

As a consequence, on short timescales the extent of the cusp built up by the IMBH shall follow the the evolution of the core radius, thus expanding or shrinking according to energy injection to the half-mass radius expansion due to few-body interactions. On longer timescales, the size of the IMBH sphere of influence would tend to roughly the same value both for clusters with concentrated and shallow profiles, following the behaviour of $r_{c}$. As discussed in the next section, this would also affect the rate of high-velocity ejections.

3.2 High velocity stars

A binary star of total mass $m=m_{1}+m_{2}$ and semi-major axis $a$ is tidally separated whenever it approaches the IMBH within the tidal radius (Hills, 1988; Yu & Tremaine, 2003)

[TABLE]

The typical (average) velocity of the ejected star is (Bromley et al., 2006)

[TABLE]

For an unequal-mass binary, due to momentum conservation the ejection speeds of the primary and secondary are

[TABLE]

respectively. While weakly dependent on the IMBH mass, the ejection velocity depends mostly on the binary semi-major axis: stars that were bound in tighter binaries are ejected with larger velocities. For what concerns the ejection rates, binaries can be disrupted in the full loss-cone or the empty loss-cone regime, according to the mass of the IMBH and the velocity dispersion of the surrounding environment (Merritt, 2013). In the full loss-cone regime, binary stars are scattered in and out of the loss-cone along their orbital path, while any star deflected into the loss-cone is disrupted within a dynamical time, in the empty loss-cone regime. In the case of binary disruptions by an IMBH in the core of a star cluster, Pfahl (2005) found a typical rate for the full loss-cone regime

[TABLE]

while, in case of empty loss-cone

[TABLE]

Here, $\sigma$ is the core velocity dispersion and $n$ the number density of the host cluster, and $f_{b}$ is the binary fraction. Note that while the rate depends on the binary semi-major axis in the full loss-cone regime, it is independent of in the empty loss-cone. For typical clusters and IMBH masses, the rates are $\sim 0.1$ – $1f_{b}\ \mathrm{Myr}^{-1}$.

Our numerical integrations show that only in the model BH-BIN, i.e. when we include both primordial binaries and the IMBH, a considerable amount of high velocity escaping stars is being produced ($\approx 45$ stars with $v_{\mathrm{esc}}\geq 50\,\,\rm km\,s^{-1}$ in $800\,{\,\rm Myr}$). Other models give on average less than one star with $v_{\mathrm{esc}}\geq 50\,\,\rm km\,s^{-1}$. This fact is a strong indication for that the primary mechanism that accelerates stars to high velocities is a separation of a binary star in the tidal field of the IMBH, as previously discussed.

Fig. 2 shows distribution of the escaping stars in model BH-BIN in the time and velocity domain. Left panel displays cumulative distribution in time, i.e. it effectively provides escape rates per unit time. The escape rate is $\approx 2.5$ stars per $100\,{\,\rm Myr}$ for $T\lesssim 200\,{\,\rm Myr}$. Then we observe growth of the escape rate which is likely related to the buildup of the stellar cusp around the IMBH. After $T\approx 300\,{\,\rm Myr}$ it is nearly constant with a value of $\approx 7$ stars per $100\,{\,\rm Myr}$ (the escape rate is approximately four times larger if we push the lower limit of escaping stars to $25\,\,\rm km\,s^{-1}$). Roughly one half of stars with $v_{\mathrm{esc}}>50\,\,\rm km\,s^{-1}$ are ordinary stars while the second half accounts for compact objects. Considering the fact that compact objects represent only $\approx 11$ per cent of the total number of ‘stars’ in the cluster, they are being ejected more effectively than the lighter ordinary stars. This is likely due to that the compact objects sink to the cluster centre by means of dynamical friction where the accelerator, the IMBH, resides.

Velocities of escaping stars span over more than one order of magnitude, starting from $\approx 10\,\,\rm km\,s^{-1}$ while the highest velocity achieved in our models is slightly above $300\,\,\rm km\,s^{-1}$. For $v_{\mathrm{esc}}\gtrsim 30\,\,\rm km\,s^{-1}$ the distribution of velocities may be approximated with a power-law, $n(v_{\mathrm{esc}})\propto v_{\mathrm{esc}}^{-3}$. Combining the analytic approximations to the output of the numerical model we obtain an estimate of the rate of escaping stars per unit velocity and unit time for $v_{\mathrm{esc}}\gtrsim 30\,\,\rm km\,s^{-1}$ as

[TABLE]

As discussed, the typical ejection velocity depends on the initial properties of the binary population, while it weakly varies with the IMBH mass. Thus, changing the IMBH would only slightly affect the typical ejection speed of the high-velocity stars. On the other hand, the choice of the binary semi-major axis distribution is fundamental in determining the typical average ejection speed of the stars ejected from the cluster (Perets & Šubr, 2012; Šubr & Haas, 2016), being $v_{ej}\propto a^{-1/2}$ (Eq. 3). In our simulations, we adopted a log-normal distribution which peaks at a semi-major axis corresponding to an orbital period of $\sim 180$ yr (Duquennoy & Mayor, 1991). An initial distribution that would favour smaller semi-major axis, as e.g. a log-uniform distribution (Öpik’s law), would lead to larger ejection velocities. While not affecting the ejection velocity, the mass of the IMBH sets the typical disruption rate of binaries and, as a consequence, of high-velocity escapers, as seen in Eqs. (5) – (6). In our models, we found a total ejection rate of the order of $\sim$ one star (or compact stellar remnant) per $\sim 10\,{\,\rm Myr}$, roughly consistent with the previous equation. Fixed the cluster properties, the rate is $\propto M_{\mathrm{IMBH}}^{4/3}$ and $\propto M_{\mathrm{IMBH}}^{3}$ in the full and empty loss-cone regime, respectively. A larger IMBH mass would give larger rates, thus enhancing the number of high-velocity stars ejected from the cluster. Similarly, clusters with larger central density and smaller velocity dispersion would produce a larger amount of ejected stars. This is in turn related to the characteristics of the cluster density profile, upon which a cusp is built by the IMBH, as discussed in the previous section (see also Trenti et al., 2007). We stress that the cluster binary fraction plays a crucial role in determining the rate of ejected stars, since the rate is lineraly dependent on it. Cluster harbouring an IMBH but with a small binary fraction would give a few ejections over the cluster lifetime, as found in our model BH. Finally, we note that $f_{b}$ is not constant throughout the cluster evolution, but it decreases in time. Trenti et al. (2007) found that typically the binary fraction decreases to $\sim 60$-$80\%$ of the initial value, both for shallow Plummer and concentrated King profile. As a consequence, the ejection rate of stars is expected to decrease with time, following the evolution of $f_{b}$.

3.2.1 Population of ejected stars in host galaxy

Although we have studied a limited number of clusters, we use the inferred rate (Eq. 7) as a proxy to estimate the statistical properties of the population of stars ejected from GCs in a Milky Way-like galaxy. Even though the most important features of ejected stars are captured by Eqs. (3) – (6), further studies that include a wider span of initial conditions, importantly the cluster and IMBH size and the distribution function of binary semi-major axis, are necessary to study in detail how different clusters populate a Milky Way-like (or another) galaxy with high-velocity stars.

First, we investigate ejection of high-velocity stars from a single cluster hosting an IMBH that moves in the galactic gravitational field. We characterize the cluster orbit by means of three parameters (Fragione & Gualandris, 2018a): (i) semi-major axis $a_{\rm GC}$, (ii) eccentricity $e_{\rm GC}$ and (iii) relative inclination of the cluster orbital plane and of the Galactic disc $\eta_{\rm cl}$ ($\eta_{\rm cl}=0^{\circ}$ corresponds to an orbital plane coinciding with the Galactic disk).

The Milky Way-like galaxy potential is described with a 4-component model $\Phi(r)=\Phi_{\rm BH}+\Phi_{\rm b}(r)+\Phi_{\rm d}(r)+\Phi_{\rm h}(r)$ (Kenyon et al., 2014; Fragione & Loeb, 2017), where:

• •

$\Phi_{\rm BH}$ is the contribution of the central SMBH,

[TABLE]

with mass $M_{\rm BH}=4\times 10^{6}\ {\rm M}_{\odot}$;

• •

$\Phi_{\rm b}$ is the contribution of the spherical bulge,

[TABLE]

with mass $M_{\rm b}=3.76\times 10^{9}\ {\rm M}_{\odot}$ and scale radius $a=0.10{\,\rm kpc}$;

• •

$\Phi_{\rm d}$ accounts for the axisymmetric disc,

[TABLE]

with mass $M_{\rm disk}=5.36\times 10^{10}\ {\rm M}_{\odot}$, length scale $b=2.75{\,\rm kpc}$ and scale height $c=0.30{\,\rm kpc}$;

• •

$\Phi_{\rm halo}$ is the contribution of the dark matter halo

[TABLE]

with $M_{\rm DM}=10^{12}\ {\rm M}_{\odot}$ and length scale $r_{\rm s}=20{\,\rm kpc}$.

The potential parameters are set so that the Galactic circular velocity is $235\ {\rm km\,s}^{-1}$ at the Sun’s distance ($8.15{\,\rm kpc}$).

To generate mock populations of ejected stars, we follow Kenyon et al. (2014) prescriptions. Moreover, we assume that the IMBH ejects stars at a constant rate along the cluster orbit. When we eject a star, we randomly draw an ejection time $t_{\rm ej}$ and an observation time $t_{\rm obs}$ between zero and the star’s main-sequence lifetime. In this approach, we both account for the fact that the star interacts with the IMBH and is observed before it evolves off the main-sequence. If $t_{\rm ej}<t_{\rm obs}$, we assign the star an ejection velocity $v_{\rm esc}$ sampled from equation (7). Since we are interested in how the ejected stars populate the host galaxy, we combine the star ejection velocity with the cluster orbital velocity at the moment of ejection. We then integrate the star orbit across the host galaxy up to a maximum time $\mathcal{T}=t_{\rm obs}-t_{\rm ej}$. If at any time the star passes the virial radius (250 kpc), we consider the star ejected from the galaxy and we remove it from our calculation.

Fig. 3 illustrates the final distribution of velocity (left) and position (right) in Galactic rest frame for stars ejected from a cluster with different orbital semi-major axis when $e_{\rm cl}=0$ and $\eta_{\rm cl}=0^{\circ}$ (top), eccentricities when $a_{\rm cl}=10$ kpc and $\eta_{\rm cl}=0^{\circ}$ (centre) and $\eta_{\rm cl}$ when $a_{\rm cl}=10$ kpc and $e_{\rm cl}=0$ (bottom). The relative inclination $\eta_{\rm cl}$ between the cluster orbit and the Galactic disc does not play an important role in determining the final shape of the velocity and position distribution. On the other hand, the cluster orbital semi-major axis affects the final spatial distribution, which is peaked around the cluster semi-major axis, but not the final velocity distribution that results peaked near the cluster orbital velocity. Finally, the cluster eccentricity affects both of the distributions in velocity and position: the larger the eccentricity the larger the broadening around the peaks determined by the cluster orbital semi-major axis.

To understand how the stars ejected from a distribution of clusters populate the host galaxy, we use the previous scheme by considering a population of 100 clusters with different semi-major axis, eccentricities and inclinations. We sample cluster semi-major axis from a log-uniform distribution (between $a_{\rm min}=5{\,\rm kpc}$ and $a_{\rm max}=50{\,\rm kpc}$), eccentricities from a thermal distribution and inclinations from an isotropic distribution (uniform in $\cos\eta_{\rm cl}$). Fig. 4 illustrates the final velocity (left) and spatial (right) distribution of the ejected stars. Most of the ejected stars are located within $30$ kpc and have peak velocity $\approx 250\,\rm km\,s^{-1}$, with a tail extending up to $\approx 500\,\rm km\,s^{-1}$. Our model does not produce a large population of stars with extreme velocities (as a consequence of the assumed initial distribution for binary periods). Hence, we can not account for most of the unbound HVS population observed in the Milky Way (Brown et al., 2014; Brown, 2015). However, we note that the main feature of our model is to produce fast moving stars that do not point back towards the Galactic Centre, so stars have a significant non-radial component of motion (in Galactic coordinates), as in the case of runaways and hyperunaways produced in the Galactic disc (Silva & Napiwotzki, 2011; Palladino et al., 2014; Zhong et al., 2014; Vickers et al., 2015). Interestingly, recent analysis of Gaia111http://sci.esa.int/gaia/ data showed that only a few of the confirmed and candidate HVSs have orbits that can be traced back to the Galactic Centre (Boubert et al., 2018; Brown et al., 2018; Marchetti et al., 2018).

Studying RS and HVS kinematic and spectroscopic data can help constraining their origin and possibly the presence of an IMBH. Spectroscopic data can be used to get radial velocities, while tangential velocities can be obtained from proper motion measurements. The combination of radial and tangential velocity, along with the position in the sky of the stars, gives the full 6-D phase space information to study the stars’ orbits. Hoogerwerf et al. (2001) used milli-arcsecond accuracy astrometry from Hipparcos and from radio observations of the orbits of $56$ RSs and nine compact objects with distances $\lesssim 700$ pc, to identify their parent stellar group. Heber et al. (2008) studied the mass, evolutionary lifetime and kinematics of HD 271791 by using proper motion measurements from a collection of catalogues. They found that the likely birthplace is the outer Galactic disc, while the Galactic Centre is ruled out. Recently, Lennon et al. (2019) and Renzo et al. (2019) used Gaia data to show that the dynamics of the very-massive runaways VFTS 16 and VFTS682 are consistent with them having been ejected from the young massive cluster R136. Finally, also Hattori et al. (2018) used Gaia data to study the origin of hyper-runaway subgiant LAMOST-HVS1, and found that it was likely ejected from near the Norma spiral arm dynamically as a consequence of a few-body encounter or a Hills ejection by an IMBH. Thanks to the high precision of Gaia proper motion similar studies can be performed for known and candidate high-velocity stars, that might have been ejected through the mechanism discussed in this work, thus disclosing new IMBH candidates.

3.3 Low velocity escapers

While the stars ejected from the star clusters at high velocities may have a specific appeal, their number counts are relatively small and, therefore, it may be difficult to make direct relations of high velocity stars and star clusters in observational data. In this situation, their low velocity counterparts may serve as an important diagnostic tool as they are more numerous and are found closer to their parent clusters. Our simulations indicate (see Fig. 5) that star clusters hosting an IMBH and a considerable population of primordial binaries produce larger number of stars escaping from the cluster with velocities $1\,\,\rm km\,s^{-1}\lesssim v_{\mathrm{esc}}\lesssim 30\,\,\rm km\,s^{-1}$ than their counterparts without an IMBH or a primordial binary population. In Fig. 6 we plot the temporal evolution of the number of stars with velocities greater than $1.5v_{\mathrm{esc}}$ as a rough measure of the cluster evaporation rate. Here, $v_{\mathrm{esc}}$ is determined by a method developed by Hénon (1971) for Monte Carlo simulations of spherically symmetric star clusters: At the radius $r_{i}\leq r\leq r_{i+1}$, with $r_{i}$ being ordered radial distances of stars from the cluster centre, the escape velocity is determined from the mean potential,

[TABLE]

where $G$ stands for the gravitational constant and $m_{i}$ are the masses of the individual stars. Fig. 5 indicates that this method approximates the mean potential of not exactly spherical $N$-body system quite well, although it cannot determine exactly whether a particular star is bound to the star cluster or not. Hence, the factor of $1.5$ used for the detection of evaporating stars. Yet another source of possible mis-detections are binary systems – for those being regularised in the NBODY6 code, we consider their centre of mass velocities. For weakly bound binaries, however, their full orbital velocities are used, which may be higher than $1.5v_{\mathrm{esc}}$ although the binary system itself may still be bound to the star cluster.

Assuming that the primordial binary population is present in real star clusters, more or less rich tidal tails of GCs may serve as an indicator for presence of the IMBH.

3.4 Direct interactions with the IMBH

In our numerical setup, we treat ‘stellar’ collisions in a quite simplified way, considering only two different classes: collisions that involve the IMBH and collisions that do not. For the collisions that do not involve the IMBH, the physical radii of the stellar mass objects are initially set up as $R_{\ast}=0.5\,\rm R_{\odot}(M_{\star}/\,\rm M_{\odot})$. If the radial separation of any of two stellar mass objects becomes less than the sum of their radii, they are merged into one object whose mass is the sum of masses of the merging objects (i.e. without any mass loss). Up to the factor $0.5$, the relation used here corresponds to a classical formula for radii of low mass stars (see e.g. Lang, 1974). The factor $0.5$ in the present prescription decreases the distance that leads to a merger, and is to ensure that the encountering stars really merge instead of merely having a close passage. Note that while this prescription may be accepted as a raw approximation for collisions of main sequence stars, it is definitely not adequate for stellar mass compact objects (WDs, NSs, BHs). Still, this is not a big issue, as even though there are collisions among compact objects occurring in our calculations, they are not frequent enough to considerably alter the mass spectrum and, consequently, the overall results.

For the collisions that do involve the IMBH, the criterion that distinguishes a merger from a fly-by is that the stellar mass object (i.e. a star or a stellar remnant) comes closer to the IMBH than the tidal radius of a main-sequence star of that mass. In that case, similar to the case of encounters without the IMBH, it is swallowed completely by the IMBH, regardless of its type and without any mass loss. Again, this prescription may be considered as very raw, but still plausible for main sequence stars (see, e.g., Guillochon & Ramirez-Ruiz, 2013, for more rigorous description of tidal disruption events). However, it may not be adequate for compact objects, as the real tidal disruption radii for both WDs and NSs are orders of magnitudes smaller than the tidal radii of main-sequence stars, which is considered as the decisive criterion here. For NSs and stellar mass BHs, one would also have to consider gravitational radiation which would help to drag them towards the IMBH, but post-Newtonian dynamics is switched off in our calculations. Despite this caveat, we plot growth of the IMBH through mergers with stars and compact objects for models BH and BH-BIN in left panel of Fig. 7. What may be interesting and worth studying with more elaborated models is the fact that the rate of growth of the IMBH roughly correlates with ejection rate of high velocity stars. Not only it saturates at a constant rate which at least by the order of magnitude corresponds to the ejections in model BH-BIN, but the correlation also goes across the models in that model BH-BIN not only produces the most high velocity stars, but also exhibits the considerable growth of the IMBH, while essentially no growth is found for model BH.

Thus, in our simple setup, collision rates of stars / compact objects with the IMBH are rather straightforwardly related to the growth rate of the IMBH. The right panel of Fig. 7 shows the cumulative distribution of close encounters with the IMBH for model BH-BIN distinguished for stars and compact objects. Trends found in this plot qualitatively agree with those found in the left panel of Fig. 2, i.e., the ejection rate as well as the rate of close encounters is initially smaller for the compact objects than for the main-sequence stars, while their relation is opposite at the end of the integrations. In the case of star–IMBH interactions, the close encounter rate may be interpreted as a rate of potentially observable tidal disruption events. The interpretation is much less straightforward for compact objects as those would require some mechanism that brings them even closer to the IMBH (e.g., two-body relaxation and/or emission of gravitational waves) in order for them to merge with the IMBH. For example, white dwarfs that dominate population of compact stellar-mass objects in real star clusters have tidal radii approximately by a factor of $100$ smaller than what was considered in our simulations.

However, all models that we consider imply that the IMBH grows only slowly through stellar dynamical processes in the environment that we simulate (a moderately old and moderately massive star cluster). This means that an IMBH is unlikely to grow from a small seed through stellar dynamics in this setting. This does not exclude that other conditions may be more favourable for the rapid growth of an IMBH, like for instance in very young clusters that still contain massive stars. Massive stars are more extended and therefore have a larger cross section for collisions. Portegies Zwart & McMillan (2002) found indeed in their $N$-body simulations of young star clusters that run-away mergers of massive stars can lead to objects that qualify as IMBH-progenitors. It may also be possible that a substantial growth of the IMBH can continue until later times in much more dense and massive clusters, like massive GCs and UCDs.

4 Discussions and Conclusions

In this paper, we have discussed the dynamical interaction of a GC of moderate mass and density with an intermediate-mass black hole. The main advancement of this study in comparison to earlier studies is that we include simultaneously to the IMBH a large population of binaries into our models for the GC, which makes the computations very demanding. Our main results are:

• •

All models presented in this paper started with identical density and velocity dispersion profiles given by the Plummer distribution. Clusters hosting the IMBH evolve towards power-law density profile with slope $\approx-7/4$ which is practically identical for both cases (i.e. regardless of the stellar multiplicity). At the end of our integrations ($T\approx 0.8\,{\,\rm Gyr}$) the profiles of cluster with and without IMBH considerably differ for $r\lesssim r_{\mathrm{h}}$; most prominent is the difference below $0.1\,r_{\mathrm{h}}$.

• •

In our calculations, close interactions with the IMBH led to direct increase of its mass (and removal of the interacting star from the integration). This approach definitely overestimates the growth rate of the IMBH, i.e. our results can only serve as an upper estimate for evolution of mass of the IMBH. Considering that in the most prominent case, model BH-BIN with the large binary fraction, the IMBH grew from $1000\,\,\rm M_{\odot}$ to $\approx 1050\,\,\rm M_{\odot}$ over $0.8\,{\,\rm Gyr}$. This indicates that at least for lightweight GCs, the collisions of stars with the IMBH are not able to increase its mass by a factor of two on cosmological time-scales.

• •

The most prominent feature of the model which included both the IMBH and large fraction of stellar binaries is ejection of high-velocity stars. As the ejection rate was about two orders of magnitude smaller for models either with the IMBH but no primordial binaries or with primordial binaries but no IMBH, it is clear that the acceleration is due to interaction of binary stars with the IMBH through the Hills (1988) mechanism. The maximum velocity we observed in our model was slightly above $300\,\,\rm km\,s^{-1}$, but this was rather exceptional case. The distribution of velocities of stars escaping from the cluster could be approximated by a power-law with a slope of $\approx-3$ (see equation 7) for $v_{\mathrm{esc}}\gtrsim 30\,\,\rm km\,s^{-1}$. Total ejection rate in this velocity interval is of the order of one star (or compact stellar remnant) per $10\,{\,\rm Myr}$. This rate, together with the fact that these stars pass a distance $>300\,{\,\rm pc}$ over $10\,{\,\rm Myr}$, indicates that there is quite a low probability of finding the high velocity escapers in the vicinity of parent cluster, i.e. serving possibly as indicators of the presence of the IMBH.

• •

The high velocity stars ejected from star clusters hosting an IMBH would contribute to the overall distribution of high-velocity stars in the host galaxy. The composite distribution of high-velocity stars from a population of cluster orbiting the host galaxy (semi-major axis from a log-uniform distribution between $a_{min}=5$ kpc and $a_{max}=50$ kpc and thermal eccentricities) generate stars that are mostly within $\sim 30$ kpc with peak velocity of $\sim 250\,\,\rm km\,s^{-1}$ and a tail extending up to $\sim 500\,\rm km\,s^{-1}$. Since our model does not produce a large population of stars with extreme velocities, we can not account for most of the unbound HVS population observed in our Galaxy, but our mechanism can produce runaway and hyperunaway stars.

• •

The models of star clusters that differ just by presence of the IMBH (BIN vs BH-BIN) considerably differ by amount of stars being accelerated above the escape velocity (see Figs. 5 and 6). Beside the high velocity escapers which are relatively rare, model with the IMBH (BH-BIN) produces a relatively numerous population of low velocity escapers with velocities $\gtrsim 1\,\,\rm km\,s^{-1}$. These form a diluted extended halo of the cluster which may be potentially observable for sufficiently isolated star clusters (such that cluster stripping by the galactic tidal field does not overlay the evaporation caused by the cluster internal dynamics).

Our current results are restricted to star clusters of moderate mass and density and also cover only a fraction of their life-time. This limitation stems from a large numerical complexity of our models and moderate computational resources available (one integration of the model BH-BIN took approximately half a year of computer time on a 40 CPU core sytem). From this point of view, our results have to be taken as another step on the way towards fully realistic numerical models of GCs hosting an IMBH. Beside predicting properties of such systems, our results may also serve for calibration of other approaches, e.g. semi-analytic calculations (Fragione, Ginsburg & Kocsis, 2018; Fragione, Leigh, Ginsburg & Kocsis, 2018), scattering of individual binaries on the IMBH in isolated three-body approximation (Fragione & Gualandris, 2018a) or Fokker-Planck methods of integration of star clusters. We finally stress that our modelisation of the high-velocity stars ejected from GCs is only a proxy to the real population. As discussed, the main properties of the ejected stars may be inferred from Eqs. (3) – (6), but further studies that consider different cluster and IMBH sizes and binary properties are highly desirable to precisely model the ejected population of stars and compare it to other scenarios (Fragione & Gualandris, 2018a).

The mass of the cluster considered in our simulations together with an assumption of the canonical IMF implies less than 100 stellar remnants to be stellar-mass black holes. Despite the mass function of newly born stellar black holes is unknown as it is a subject to many processes (the stellar collapse itself, subsequent growth through mergers, ejection due to natal kicks and mutual scattering), we may expect that there may be no more than 10 black holes of mass $\gtrsim 10\,\,\rm M_{\odot}$. This number is small enough to assume that this component of the cluster will not affect considerably its overall distribution. Nevertheless, it is likely that these massive stellar black holes will sink towards the IMBH and will dominate the central cusp (which is formed by no more than several tens of stars in our models BH-BIN and BH). Hence, it is worth to discuss whether these massive black holes may influence the processes that take place in the vicinity of the IMBH, namely the Hill’s process. In order to shed some light onto this problem, we have evaluated properties of orbits of binaries before one of their components was accelerated to high velocity. It appears that semi-major axes of these systems lie mostly in the range $\sim 0.1\,{\,\rm pc}$-$1\,{\,\rm pc}$, i.e. they are just penetrating into the cusp on radial orbits with high velocities and they spend only short time there. It is then likely that there will be a rather small probability of strong interaction with massive black holes and, therefore, we suggest that our results regarding ejection of stars from the cluster via interactions with the IMBH will not be considerably affected. What may be affected is observability of the central cusp which, if dominated by black holes, will not be detectable by conventional observations. (Note, however, that the cusp formed by several tens of stars will be definitely hard to observe anyway.)

We did not included an external tidal field, but we do not expect any significant changes regarding the processes in the vicinity of the IMBH, i.e. the formation of the cusp and the production of high velocity escapers. Only the presence of a tidal field would increase the evaporation of the cluster and the escape rate of stars with small velocities ($\lesssim 30\,\rm km\,s^{-1}$).

5 Acknowledgements

We thank Sverre J. Aarseth for helpful discussions on optimizations of the NBODY6 code. This work was supported by the Czech Science Foundation through the project of Excellence No. 14-37086G. Most of the numerical calculations were performed on the computational cluster Tiger at the Astronomical Institute of the Charles University. GF acknowledges hospitality from Ladislav Šubr and Charles University of Prague, where the early plan for this work was conceived. GF is supported by the Foreign Postdoctoral Fellowship Program of the Israel Academy of Sciences and Humanities. GF also acknowledges support from an Arskin postdoctoral fellowship and Lady Davis Fellowship Trust at the Hebrew University of Jerusalem.