Intermediate-mass black holes in binary-rich star clusters Free

Abstract

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 \mathit{ 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, 2017; Fragione & Kocsis 2018). However, for the intermediate-size population, the so-called intermediate-mass black holes (IMBHs) (⁠|$100\, \rm M_\odot \lesssim \mathit{ 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 modelling, namely ω Cen (Noyola et al. 2010; Jalali et al. 2012; Baumgardt 2017) and 47 Tuc (Kızıltan, Baumgardt & Loeb 2017), but other authors (Anderson & van der Marel 2010; van der Marel & Anderson 2010; Zocchi, Gieles & Hénault-Brunet 2017, 2019) 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 ∼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_SMBH–σ relation (σ is the velocity dispersion) also holds 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, Gürkan & Rasio 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 ∼0.1 per cent 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, Hernquist & Loeb 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 modelling 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 et al. 2018; Fragione, Ginsburg & Kocsis 2018; Arca-Sedda & Gualandris2018; Arca-Sedda & Capuzzo-Dolcetta 2019). 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 ∼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 (Amaro-Seoane & Santamaría 2010; Leigh et al. 2014; Fragione et al. 2018; Fragione & Leigh 2018a). Besides GW events, TDEs may also take place in the centre of a cluster hosting an IMBH (Baumgardt, Makino & Ebisuzaki 2004a,b; Fragione et al. 2018), whose rate (10⁻⁶–10⁻³ yr⁻¹) can be of the same order of the TDEs driven by the host galaxy SMBH.

Several efforts have been made in modelling GCs with IMBHs through direct N-body simulations. Baumgardt, Makino & Hut (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 per cent 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, Baumgardt & Kruijssen (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_IMBH/M_GC (where M_IMBH and M_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, Trenti & Hut 2006; Trenti, Heggie & Hut 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, Kobayashi & Rossi 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 per cent). 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 to 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; Sesana, Haardt & Madau 2006; Gualandris & Portegies Zwart 2007; Zubovas, Wynn & Gualandris 2013; Capuzzo-Dolcetta & Fragione 2015; Fragione & Capuzzo-Dolcetta 2016; Šubr & Haas 2016; Boubert et al. 2017; Fragione, Capuzzo-Dolcetta & Kroupa 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 50 000 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, 50 000 particles sounds like a modest number in comparison to the 10⁶ particles considered by Wang et al. (2016), but we note that in the simulations of Wang et al. (2016) only 5 per cent of the stars are in primordial binaries (in contrast to 50 per cent 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 set-up for the system that is at the centre of attention in this 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 α = −1.3 for |$M_\star \lt 0.5\, \rm M_\odot$| and α = −2.3 for |$0.5 \, \rm M_\odot \le M_\star \lt 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 Myr at most, so that the birth of all stars in an initial starburst 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 50 000. Of these particles, 44 500 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 clusters.

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 44 500, and the upper mass limit of the IMF, m_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\, \,$| 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 figs 12 and 16 in Woosley, Heger & Weaver 2002), but their numbers are suppressed strongly by the steep slope of the canonical IMF at high masses.

We also finally include 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 are 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 this 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_bin, is defined as the number of binaries over the number of centre-of-mass systems, i.e. binaries and single stars in our case. We set f_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_bin = 1, which was suggested by Kroupa 1995 for dynamically unevolved star clusters). We note that Wang et al. (2016) assumed f_bin = 0.05 in their N-Body simulation with 10⁶ individual stars, while Sollima et al. (2007) found f_bin > 0.06 for the low-density GCs they studied, and in some cases f_bin ≈ 0.5 seems likely. This motivates our high choice for f_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_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 semimajor axes of primordial binaries follow the distribution found by Duquennoy & Mayor (1991) for the Galactic field in all models we consider, i.e. a lognormal distribution which peaks at a semimajor axis corresponding to an orbital period of ∼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 set-up 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 50 000 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, Kroupa & Baumgardt 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_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 are given in Table 1. For each model, seven realizations 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 realizations 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 neighbour 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 neighbour lists of all other (star) particles. The modification of the neighbour 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 \ge 100{\, \rm Myr}$| (see Fig. 1). The build-up 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 realizations of a particular model may differ from each other. Once the central cusp approaches the equilibriums state, its evolution slows down and both models which 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 ϱ(r) ∝ 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).

Figure 1.

Radial density profiles of models BH-BIN (red crosses), BIN (green circles), and BH (blue triangles) at six different epochs (indicated in each panel); dotted line indicates the power law with slope index −7/4.

Open in new tabDownload slide

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 semimajor axis smaller than |${\lesssim } 10^{-4}{\, \rm pc}$| and eccentricity e > 0.99 (see also MacLeod, Trenti & Ramirez-Ruiz 2016, for discussion of dynamical formation and properties of close stellar companions to an IMBH). Peak velocities of these stars during passages through the per centres of their orbits are few thousands of kilometres 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 regularization 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 per cent 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.

3.2 High-velocity stars

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 (≈45 stars with |$v_\mathrm{esc}\ge 50\, \, \rm km\, s^{-1}$| in |$800\, {\, \rm Myr}$|⁠). Other models give on average less than one star with |$v_\mathrm{esc}\ge 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-hand panel displays cumulative distribution in time, i.e. it effectively provides escape rates per unit time. The escape rate is ≈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 build-up of the stellar cusp around the IMBH. After |$T \approx 300\, {\, \rm Myr}$| it is nearly constant with a value of ≈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}\gt 50\, \, \rm km\, s^{-1}$| are ordinary stars while the second half accounts for compact objects. Considering the fact that compact objects represent only ≈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.

Figure 2.

Properties of ejected stars from model BH-BIN. Left: cumulative number of stars escaping with velocities |$\gt 50\, \, \rm km\, s^{-1}$|⁠. Green circles correspond to main-sequence stars while blue triangles show escape rate of compact objects. Red crosses represent the total and the thin dotted line indicates escape rate of 7 stars per |$100\, {\, \rm Myr}$|⁠. Right: number counts of escaping stars up to |$T = 800\, {\, \rm Myr}$| in logarithmically equal-width bins, i.e. the displayed quantity is proportional to |$v_\mathrm{esc}\, n(v_\mathrm{esc})$| with n(⁠|$v$|_esc) being the distribution function of escape velocities. For |$v_\mathrm{esc}\gtrsim 30\, \, \rm km\, s^{-1}$| the displayed quantity may be approximated with a power law |$\propto v_\mathrm{esc}^{-2}$| (indicated by thin dotted line).

Open in new tabDownload slide

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 semimajor 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 ∝ a^−1/2 (equation 3). In our simulations, we adopted a lognormal distribution which peaks at a semimajor axis corresponding to an orbital period of ∼180 yr (Duquennoy & Mayor 1991). An initial distribution that would favour smaller semimajor 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 equations (5)–(6). In our models, we found a total ejection rate of the order of ∼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 linearly 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 ∼60–|$80{{\ \rm per\ cent}}$| 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 (equation 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 equations (3)–(6), further studies that include a wider span of initial conditions, importantly the cluster and IMBH size and the distribution function of binary semimajor 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) semimajor axis a_GC, (ii) eccentricity e_GC, and (iii) relative inclination of the cluster orbital plane and of the Galactic disc η_cl (η_cl = 0° corresponds to an orbital plane coinciding with the Galactic disc).

The Milky Way-like galaxy potential is described with a four-component model Φ(r) = Φ_BH + Φ_b(r) + Φ_d(r) + Φ_h(r) (Kenyon et al. 2014; Fragione & Loeb 2017), where:

• Φ_BH is the contribution of the central SMBH,

  \begin{eqnarray*} \Phi _{\rm BH}(r)=-\frac{GM_{\rm BH}}{r}, \end{eqnarray*}

  (8)

  with mass M_BH = 4 × 10⁶ M_⊙;
• Φ_b is the contribution of the spherical bulge,

  \begin{eqnarray*} \Phi _{\rm b}(r)=-\frac{GM_{\rm bul}}{r+a}, \end{eqnarray*}

  (9)

  with mass M_b = 3.76 × 10⁹ M_⊙ and scale radius |$a=0.10{\, \rm kpc}$|⁠;
• Φ_d accounts for the axisymmetric disc,

  \begin{eqnarray*} \Phi _{\rm disc}(R,z)=-\frac{GM_{\rm disc}}{\sqrt{(}R^2+(b+\sqrt{c^2+z^2})^2)}, \end{eqnarray*}

  (10)

  with mass M_disc = 5.36 × 10¹⁰ M_⊙, length scale |$b=2.75{\, \rm kpc}$| and scale height |$c=0.30{\, \rm kpc}$|⁠;
• Φ_halo is the contribution of the dark matter halo

  \begin{eqnarray*} \Phi _{\rm halo}(r)=-\frac{GM_{\rm DM}\ln (1+r/r_s)}{r}. \end{eqnarray*}

  (11)

  with M_DM = 10¹² M_⊙ 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_ej and an observation time t_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_ej < t_obs, we assign the star an ejection velocity |$v$|_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 semimajor axis when e_cl = 0 and η_cl = 0° (top), eccentricities when a_cl = 10 kpc and η_cl = 0° (centre), and η_cl when a_cl = 10 kpc and e_cl = 0 (bottom). The relative inclination η_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 semimajor axis affects the final spatial distribution, which is peaked around the cluster semimajor 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 semimajor axis.

Figure 3.

Arbitrarily normalized final distributions of velocity (left) and position (right) in Galactic rest frame for stars ejected from a cluster of with different orbital semimajor axis when e_cl = 0 and η_cl = 0° (top), eccentricities when |$a_{\rm cl}=10{\, \rm kpc}$| and η_cl = 0° (centre), and η_cl when |$a_{\rm cl}=10{\, \rm kpc}$| and e_cl = 0 (bottom).

Open in new tabDownload slide

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 semimajor axis, eccentricities, and inclinations. We sample cluster semimajor 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 η_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 cannot account for most of the unbound HVS population observed in the Milky Way (Brown, Geller & Kenyon 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, Smith & Grebel 2015). Interestingly, recent analysis of 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, Rossi & Brown 2018).

Figure 4.

Final velocity (left) and spatial (right) distribution in the Galactic rest frame for stars ejected by 100 clusters in a Milky Way like galaxy. Cluster semimajor axis is drawn from a log-uniform distribution (⁠|$a_{\mathrm{ min}}=5{\, \rm kpc}\le a_{\mathrm{ cl}}\le a_{\mathrm{ max}}=50{\, \rm kpc}$|⁠), eccentricities from a thermal distribution and inclinations from an isotropic distribution (uniform in cos η_cl).

Open in new tabDownload slide

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 6D phase-space information to study the stars’ orbits. Hoogerwerf, de Bruijne & de Zeeuw (2001) used milliarcsec accuracy astrometry from Hipparcos and from radio observations of the orbits of 56 RSs and nine compact objects with distances ≲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. (2018) 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

Figure 5.

Top panels: Velocity versus radial distance from the cluster centre at T = 800 Myr for all three models of star clusters (one randomly selected run per model). Thick green line indicates escape velocity at given radius determined by the method of Hénon (1971); thin dotted lines in each panel show escape velocity limits for all other models for comparison (the lines for models BH and BH-BIN lie on the top of each other). Bottom: number counts of stars with velocity |$v$| > 1.5|$v$|_esc in logarithmically equal-width radial bins.

Open in new tabDownload slide

Figure 6.

Temporal evolution of the number of escaping (unbound; |$v$| > 1.5|$v$|_esc) stars for models BH-BIN (solid red line), BH (dotted blue), and BIN (dashed green).

Open in new tabDownload slide

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 set-up, 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 this 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 TDEs). 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-hand 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.

Figure 7.

Left: Temporal evolution of mass of the IMBH for models BH-BIN (solid red line) and BH (dotted blue) under assumption that all close encounters of stars and compact objects contribute to the mass of the IMBH. Right: Cumulative distribution of close encounters of stars (dotted blue line) and compact objects (solid red) with the IMBH in model BH-BIN.

Open in new tabDownload slide

Thus, in our simple set-up, collision rates of stars/compact objects with the IMBH are rather straightforwardly related to the growth rate of the IMBH. The right-hand 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-hand 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 TDEs. 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 GWs) 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 slowly only 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 IMBH. 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 ≈−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 ≲ r_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\, \,$| 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 ≈−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 |$\gt 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 (semimajor 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 ∼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 cannot 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 versus 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 lifetime. 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 system). 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 et al. 2018, 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 modelization 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 equations (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 on to this problem, we have evaluated properties of orbits of binaries before one of their components was accelerated to high velocity. It appears that semimajor axes of these systems lie mostly in the range |${\sim } 0.1\,$|–|$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 include 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}$|⁠).

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.

Footnotes

REFERENCES