MOCCA Survey Database I: Dissolution of tidally filling star clusters harbouring BH subsystems

TL;DR

This paper explores a third mechanism for star cluster dissolution driven by stellar-mass black hole subsystems, leading to abrupt cluster disintegration and potential observational signatures like dark clusters and free-floating black holes.

Contribution

It introduces a new dissolution mechanism involving black hole subsystems causing rapid cluster disintegration, expanding understanding beyond traditional mass-loss processes.

Findings

Black hole subsystems can cause abrupt cluster dissolution.

Dissolution results from energy generation and tidal stripping.

Potential observational signatures include dark clusters and free black holes.

Abstract

We investigate the dissolution process for dynamically evolving star clusters embedded in an external tidal field by exploring the MOCCA Survey Database I, with focus on the presence and evolution of a stellar-mass black hole subsystem. We argue that the presence of a black hole subsystem can lead to the dissolution of tidally filling star clusters and this can be regarded as a third type of cluster dissolution mechanism (in addition to well-known mechanisms connected with strong mass loss due to stellar evolution and mass loss connected with the relaxation process). This third process is characterized by abrupt cluster dissolution connected with the loss of dynamical equilibrium. The abrupt dissolution is powered by strong energy generation from a stellar-mass black hole subsystem accompanied by tidal stripping. Additionally, we argue that such a mechanism should also work for even tidally under-filling clusters with top-heavy initial mass function. Observationally, star clusters which undergo dissolution powered by the third mechanism would look as a 'dark cluster' i.e. composed of stellar mass black holes surrounded by an expanding halo of luminous stars (Banerjee & Kroupa 2011), and they should be different from 'dark clusters' harbouring intermediate mass black holes as discussed by Askar et al. (2017a). An additional observational consequence of an operation of the third dissolution mechanism should be a larger than expected abundance of free floating black holes in the Galactic halo.

Figures

Equations

$$T_r=\frac{0.065v_m^3}{G\rho mln\Lambda },$$

$$T_{rh}=\frac{0.138M{R}_h^{3/2}}{G^{1/2}mln\Lambda },$$

Full text

††pubyear: 2018††pagerange: Mocca-Survey Database I: Dissolution of tidally filling star clusters harbouring black hole subsystems–LABEL:lastpage

1 Introduction

Understanding the evolution of star clusters, and particularly the way in which they lose mass, has been one of the long standing problems of stellar dynamics. Already in 30s and 40s of the previous century Ambartsumian (1938) and Chandrasekhar (1942) pointed out the role relaxation plays in star cluster evolution. In the 90s Chernoff & Weinberg (1990) showed that star cluster structure, stellar evolution and tidal field of the parent galaxy have a very strong influence on cluster evolution and its dissolution time. This was a seminal study, which later resulted in many more detailed papers discussing the final state of star clusters. It was shown that several other factors can also influence the lifetime of star clusters, such as: crossing time-scale (Whitehead et al., 2013), primordial binary population (e.g. Tanikawa & Fukushige, 2009), type of the Galactic orbit (e.g. Baumgardt & Makino, 2003), the form of the Galactic potential and tidal shocking (e.g. Gnedin & Ostriker, 1997), and the properties of dark remnants retained in star clusters (Banerjee & Kroupa, 2011; Contenta et al., 2015).

In this paper, inspired by the works of Fukushige & Heggie (1995); Banerjee & Kroupa (2011); Whitehead et al. (2013) and Contenta et al. (2015) where cluster dissolution is discussed from the point of view of stellar evolution mass loss and formation of a dark subcluster, we decided to take a closer look at the process of dissolution for tidally filling clusters in which a massive subsystem of stellar-mass black holes (BHS) forms and survives until the death of the cluster.

There are generally two types of cluster dissolution discussed in the literature. Slow, connected with relaxation-driven mass loss, called ’skiing’, and abrupt, considered as dynamical, connected with stellar evolution mass loss, called ’jumping’ (Contenta et al., 2015). Fukushige & Heggie (1995) attributed the final abrupt cluster dissolution to the loss of dynamical equilibrium within the cluster. In such a situation, a cluster will not undergo core collapse111In the paper, by the term ’core collapse’ we will mean the phase of the cluster evolution, in which the energy needed to support cluster structure is provided not through release of potential energy by contracting central cluster parts, but by energy generation, in the cluster core, by binaries or an intermediate mass black hole. It is worth to note that observers use the term ’core collapse’ to mean a particular spatial cluster structure that has a central power-law surface brightness profile and cannot be well fitted by a King profile, which is characteristic for ’pre-collapse’ clusters. For clusters supported by binaries formed in dynamical interactions or by an intermediate mass black hole, both definitions are equivalent. Only in the case of primordial binaries and BHS they have different meaning.. Also as it was pointed out by Banerjee & Kroupa (2011), in the case of a strong tidal field, the final stages of cluster evolution can be connected with fast removal of luminous stars and formation of a so called ’dark cluster’ consisting of BHs and neutron stars (NSs). They attributed the cluster dissolution to the loss of virial equilibrium for luminous stars and rapid dynamical evolution of a bound very low $N$ cluster composed of evolved stars.

In this paper we introduce a third mechanism for cluster dissolution. We will investigate processes responsible for the ’jumping’ cluster dissolution connected with the formation of a massive BHS and strong energy generation in dynamical interactions between BHs and BH-BH binaries (BBHs). We will examine what is the influence of the cluster initial properties on the cluster lifetime and the way in which it is dissolving.

The structure of the paper is as follows. In Section 2 we describe the used Mocca and N-body models. In Section 3 we present our results with some discussion on possible mechanisms leading to final cluster dissolution. The final Section summarizes our conclusions and tries to relate early cluster dissolution to the possible population of BHs in the Galactic halo.

2 Model

In this paper, we use the results from the Mocca-Survey Database I (Askar et al., 2017a) that contains about 2000 models of globular clusters (GC) with different initial masses and structural and orbital parameters. These models were simulated with the Mocca code (Hypki & Giersz, 2013), which treats the relaxation process using the method described by Hénon (1971), that was significantly improved by Stodolkiewicz (1982, 1986), and more recently by Giersz and his collaborators (Giersz et al., 2008; Giersz et al., 2013, 2015; Hypki & Giersz, 2013, and reference therein). Mocca code consists of the following ingredients: Sse and Bse codes (Hurley et al., 2000, 2002) for treating binary and stellar evolution, while strong binary-single and binary-binary interactions are handled by the Fwebody code (Fregeau et al., 2004), stars are escaping from the cluster according to the description given by Fukushige & Heggie (2000). Escape is no longer instantaneous, but takes place after some delay. The parameters of the cluster models in the Mocca-Survey Database I are described in Table 1 in Askar et al. (2017a). In short, for half of the simulated models, supernovae (SNe) natal kick velocities for NSs and BHs are assigned according to a Maxwellian distribution, with velocity dispersion of 265 $kms^{-1}$ (Hobbs et al., 2005). In the remaining cases, BH natal kicks were modified according to the mass fallback procedure described by Belczynski et al. (2002). Metallicities of the models were selected as follows: Z = $0.0002,0.001,0.005,0.006$ and $0.02$. All Mocca models were characterized by a Kroupa (2001) initial mass function, with a minimum and maximum initial stellar mass of $0.08M_{\odot}$ and $100M_{\odot}$, respectively. The GC models were described by King (1966) models with central concentration parameter values $W_{0}=3,6$ and $9$. They had tidal radii ($R_{t}$) equal to: $30,60$ or $120$ pc, while the ratios between $R_{t}$ and half-mass radius ($R_{h}$) were $50,25$ or the model was tidally filling. The primordial binary fractions were chosen to be: $5\%,10\%,30\%$ and $95\%$. Models characterized by an initial binary fraction equal to or lower than $30\%$ had their initial binary eccentricities selected according to a thermal distribution (Jeans, 1919), the semi-major axes according to a flat logarithmic distribution, and the mass ratio according to a flat distribution. For models containing a larger binary fraction, the initial binary properties were instead selected according to the distribution described by Kroupa (1995, 2011). The models consist of $40000$, $100000$, $400000$, $700000$ and $1200000$ objects (stars and binaries). The GCs were assumed to move on a circular orbit at Galactocentric distances between $1$ and $50$ kpc. The Galactic potential was modelled in the simple point-mass approximation, taking as central mass the value of the Galaxy mass enclosed within the GCs orbital radius. The GC rotation velocity was set to $220$ $kms^{-1}$ for the whole range of galactocentric distances. As it was pointed out by Askar et al. (2017a), the initial conditions assumed to create the Mocca-Survey Database I were not specifically selected to reproduce the Galactic GC population. Nevertheless, their observational parameters calculated at the present-day exhibit a remarkably good agreement with Milky Way GCs (see Fig.1 in Askar et al. (2017a)).

For the purpose of this paper we mainly selected tidally filling models and models with $95\%$ binary fraction. The choice of such a high value of binary fraction may seem incomprehensible to the reader, especially when the observed GC binary fractions are recalled (Milone et al., 2012, of the order of $10\%$). However, we would like to point out, that because of a very wide primordial binary period distribution (Kroupa, 1995; Belloni et al., 2017b) most of the binaries will be quickly destroyed in dynamical interactions, and final binary fraction (at the Hubble time) will be similar to the observed ones. Additionaly, as it was shown in many papers (e.g. Leigh et al., 2015; Belloni et al., 2017b, 2018, 2019) models with very high binary fractions better describe in GCs: the observed binary fractions and their spatial distributions, the observed binary properties in the field, the dependence of observed binarity as a function of binary mass, the number and spatial distribution of CVs in star clusters. Therefore, we decided to illustrate the discussion of physical mechanisms leading to abrupt cluster dissolution mostly with models with initially very high binary fraction. As it will be clearly shown in Section 3.1, also models with an initially small binary fraction are abruptly disrupted.

To check if the Mocca models properly recover the tidally filling cluster dissolution, we run an additional N-body model. The initial conditions for this model were generated from the Mocca snapshot (stars and center of binary mass: radial position, radial and tangential velocities, mass, and for binaries additionally: semi-major axis, eccentricity and second component mass) for the initial model of a cluster. The properties of the tidally filling model are: $N=100000,W_{0}=6,R_{t}=60$ pc, $R_{h}=8.8$ pc, binary fraction equal to $0.1$ and mass fallback for BH SN kicks. The model was run with the Nbody6++gpu code (Wang et al., 2015). This is a hybrid parallelized version based on the direct $N$-body code Nbody6 (Aarseth, 2003) designed for the realistic modelling of star clusters. The single and binary stellar evolution packages Sse/Bse (Hurley et al., 2000, 2002) are used222The Nbody6++gpu code version is based on the Master version (Commits on Sept 28, 2018) on GitHub: https:/1/github.com/nbodyx/Nbody6ppGPU/commits/master. The Sse/Bse packages are the original version without recent updates used in Nbody6. In this N-body code version, the SNe kicks are treated exactly as in the Mocca code.

It is important to point out that the properties of compact binaries outputted by Mocca and Nbody6++gpu strongly depend on the treatment of the common envelope phase in binaries. The version we used for our study assumes $\alpha=3$, $\lambda=0.5$ (Hurley et al., 2002), which are default values used in the Bse code. In most N-body simulations $\alpha=3$, $\lambda=0.0$ are used by default. As a result, a lot of binaries survive to the point where both components collapse into BHs, and a significant number of mergers are primordial, without any strong dynamical interactions (Belloni et al., 2017a).

3 Results

Many theoretical investigations and star cluster simulations have pointed out that for clusters embedded in a galactic tidal field the cluster dissolution depends on the initial cluster relaxation time and the concentration and shape of the cluster initial mass function (IMF). Already, Chernoff & Weinberg (1990) showed that a cluster with small concentration (King model parameter $W_{0}$ equal to 1 or 3) will very quickly dissolve (before any substantial dynamical evolution). This cluster dissolution is powered by strong stellar evolution mass loss coupled with strong cluster expansion and overflow of the $R_{t}$. Dissolution is characterized by the abrupt loss of equilibrium and very fast dissolution - ’jumping’ shape of the evolution of the total cluster mass. Models which are more concentrated are controlled by the relaxation process, not by stellar evolution mass loss. They evolve much slower and can enter the core collapse phase, provided that the relaxation was fast enough. They are characterized by the so called ’skiing’ shape of the total mass evolution. There is a clear interplay between the time-scale of the relaxation process and the time-scale of stellar evolution mass loss, which act differently depending on the cluster concentration, strength of the tidal field and IMF shape.

In this Section, we will discuss the third mechanism leading to the abrupt cluster dissolution, namely formation and strong energy generation by BHS in tidally limiting star clusters. For such models clusters undergo core collapse connected with strong mass segregation and BHS formation, and later, usually after few Gyr of balanced evolution, suddenly undergo abrupt dissolution.

3.1 Cluster dissolution - global cluster parameters

Figure 1 (top panel) shows the evolution of the fraction of cluster bound mass (normalized to the initial value) for tidally filling models with initially N=$700000$ objects, $95\%$ initial binary fraction, King model concentration $W_{0}=6$ for different tidal radii and SNe natal kicks with mass fallback ON or OFF. It is clear that models with mass fallback ON evolve much faster than models with mass fallback OFF, and are characterized with abrupt disruption, which in turn suggests loss of the virial equilibrium (Fukushige & Heggie, 1995). This behaviour is different from the case described by Contenta et al. (2015) where models with NS natal velocity kicks being zero (similar to the Mocca models with fallback ON - reduced natal kicks, but not set to zero) showed clear core collapse and dissolution controlled by the relaxation process, while models with large NS natal velocity kicks (similar to the Mocca models with fallback OFF) showed abrupt dissolution features controlled by strong mass loss due to stellar evolution. This different behaviour is connected with the maximum mass of the main sequence stars for IMF chosen by Contenta et al. (2015). Limiting the maximum mass to only 15 $M_{\odot}$ prevents formation of BHs which can segregate and form a BHS which takes over the control of global cluster evolution, as it can be seen in Section 3.2. Figure 1 (bottom panel) shows the evolution of the number of BHs bound to the system. Models with SNe natal kicks with mass fallback ON show very slow evolution of the BH number. This suggests that BHs close to the final stages of cluster evolution could play a crucial role in the determination of the cluster structure and evolution, as it can be seen in Section 3.2.

Cluster mass evolution shown in Figure 1 is very similar until the mass fraction is about 0.7 for all models irrespective of the $R_{t}$, or mass fallback option. This strongly suggests that the initial cluster evolution is governed by mass loss connected with the stellar/binary evolution of massive stars. Then the relaxation process takes over and the cluster evolution starts to depend on $R_{t}$. The larger the $R_{t}$, the longer the relaxation time and the slower the cluster evolution. From the fraction of cluster bound mass about 0.3 the evolution of cluster with mass fallback ON and OFF starts to differ. Models with fallback ON evolve faster. BHS takes over and governs the cluster evolution and finally leads to the loss of equilibrium within the cluster.

To check if the abrupt cluster dissolution shown in the Mocca simulations is real, and not connected with the physical assumptions underlying the Monte Carlo method, we decided to run a direct N-body model and check if it will show the same behaviour as the Mocca model. We have chosen a small N model with a small binary fraction, namely $N=100000$ and $10\%$ binary fraction. Such models are relatively easy for N-body simulations. The model was run with the Nbody6++gpu code (Wang et al., 2015). We used the initial Mocca snapshot (positions, velocities and masses) to setup the initial conditions for the run. As it can be seen from Figure 2, for $N=100000$ models with thick black lines, the N-body model shows the same behaviour as the Mocca model. The small differences observed from the very beginning are connected with slightly different prescriptions in the Bse code for the BH masses formed in SNe and the amount of accreted mass on BHs during collisions/mergers with other stars. The final cluster dissolution in the N-body model is slightly less abrupt than in the Mocca model, but we should remember that the Monte Carlo method is less accurate for processes which act in time-scales proportional to the dynamical time-scale. This is the case for the final cluster dissolution. Concluding, we can safely assume that the abrupt dissolution of clusters observed in the Mocca simulations are real features of cluster evolution. It is important to stress here that abrupt dissolution means a final loss of about $15\%$ of the initial cluster mass in about a few hundred Myr, compared to about dozens of Gyr of slow mass loss powered by the relaxation process.

In the literature, in order to describe the local and global cluster evolution rates, there are generally used definitions of relaxation times given in Spitzer’s book Spitzer (1987, Eqs. 2-62 and 2-63 respectively). The definition of local relaxation time is as follows:

[TABLE]

where G is the gravitational constant, $v_{m}$ is the velocity dispersion, $\rho$ is the mass density, $m$ is the average mass and $ln\Lambda$ is the Coulomb logarithm. The half-mass relaxation time if given by:

[TABLE]

where $M$ is the cluster mass. As it can be easily seen from Eqs. 1 and 2, the larger the cluster mass or the half-mass radius, the smaller the density or the average mass, the longer the relaxation time. Those dependences will be explored during discussions of the Figures presented in this Section.

In Figures 2, 3 and 4, the evolution of the fraction of cluster bound mass is shown as a function of time for models with different: number of objects, King model concentrations, binary fractions and metallicities, respectively. As it was already shown during discussion of Figure 1 most of the cluster evolution is governed by the relaxation process. Therefore, it is no surprise that models which are characterized by shorter half-mass relaxation time evolve faster and earlier enter the abrupt dissolution phase. Shorter relaxation time means: smaller number of objects (see Figure 2), smaller King concentration parameter leading to faster mass loss and stronger tidal striping (see Figure 3) and smaller binary fraction, (see Figure 4). Models with very large binary fraction, $0.95$, are initially characterized by very wide period distribution (Kroupa, 1995), so a large number of binaries are quickly disrupted in dynamical interactions and the number of objects is sharply increased leading to longer relaxation times. The dependence of the cluster disruption time on the cluster metallicity (see Figure 4) is negligible, which is not surprising. As we know from stellar evolution models, stars with larger metallicities have stronger winds and lose more mass before BH formation (Vink et al., 2001). So, BHs have slightly smaller masses and clusters should have slightly larger half-mass radii, because of enhanced stellar evolution mass loss. Larger half-mass radius means slower evolution (see Eq. 2) and smaller BH masses mean denser central parts of systems and faster removal of BHs from the system in dynamical interactions. Both effects leads to opposite behaviours, and they nearly cancel each other out. So, the abrupt cluster dissolution time does not strongly depend on metallicity.

There are only three exceptions from the picture discussed above. In Figure 2 we see that the model with the smallest number of objects, $N=40000$, does not show abrupt dissolution. This model is characterized with the shortest initial relaxation time (among models considered in this paper), about 1.6 Gyr, and the smallest number of retained BHs, about 30. In such a situation, according to Breen & Heggie (2013), BHs are quickly kicked out from the system and BHS cannot survive for a long time. So, the cluster evolution is not governed by strong energy generation by BHs in BHS. In Figure 3 we can see that the model with $W_{0}=3$ dissolves extremely fast and the model with $W_{0}=9$ does not show the abrupt dissolution feature. The very fast dissolution of tidally filling models with low King model concentration was already extensively discussed in the literature (e.g. Weinberg, 1993; Fukushige & Heggie, 1995; Whitehead et al., 2013; Contenta et al., 2015). The dissolution of such clusters is controlled by very strong initial mass loss powered by stellar/binary evolution. Relaxation process is not important at all. The situation is much different for the model with $W_{0}=9$. As it can be seen from Figure 5, the cluster enters the core collapse phase, so it has to be controlled by the relaxation process. The difference between the model with $W_{0}=6$ and the model with $W_{0}=9$ is connected with the rate of BHS evolution. Model $W_{0}=9$ is initially much denser, so its half-mass radius and half-mass relaxation time are shorter than for the $W_{0}=6$ model. Nevertheless, for both models the BH mass segregation ends nearly at the same time, about 4 Gyr. At that time, the models contain about 50 BHs and 160 BHs for $W_{0}=9$ and $W_{0}=6$, respectively. According to Breen & Heggie (2013) the evolution of BHs is controlled by the energy flow through the cluster half-mass radius, which is proportional to $E_{b}/T_{rh}\approx GM/R_{h}/T_{rh}$, where $E_{b}$ is the cluster binding energy. Models with $W_{0}=9$ have smaller half-mass radius and half-mass relaxation time than models with $W_{0}=6$, so the energy demand is larger for this model and leads to a much faster burning out of BHs - larger number of dynamical interactions leading to BH removal from the system. After the time of BH mass segregation (about 4 Gyr) the model $W_{0}=6$ contains enough BHs to form a BHS and enters the phase of balanced evolution (Breen & Heggie, 2013). Contrarily, the model with $W_{0}=9$ has too small a number of BHs and continues to remove BHs quickly to support the needed energy flow. Finally, it enters the phase when other energy sources connected with ordinary binaries take over and the cluster enters the core collapse and then core bounce phases. It is important to note that the tidal field plays an important role in the above picture. The less concentrated model, has larger $R_{h}$ and loses more mass, so it is easier for the BHS to provide the needed energy to support the cluster structure. This phase of evolution ends when the cluster suddenly loses its virial equilibrium.

There is another feature suggesting very strong influence of the BHS on the abrupt cluster dissolution. As it is shown in Figure 2, models with initially $N=700000$ and $N=1200000$ objects dissolve nearly at the same time, whereas other models with smaller $N$ follow more or less the dependence on the relaxation time - the shorter the relaxation time, the faster the cluster evolution. So, model $N=1200000$ speeds up its evolution. The faster evolution is connected with a more massive BHS and larger energy generation, which leads to stronger mass loss, additionally increased by the tidal stripping. We have to remember that the more massive cluster model has initially the same tidal and half-mass radii as the lower mass model, so the energy flow at $R_{h}$ is much larger and the BHS has to evolve faster to provide the needed energy. To a smaller extent this is also visible for the smaller N models.

At the end of this section we would like to discuss Figure 5, in which the evolution of the fraction of cluster bound mass as a function of the King concentration parameter for models with different $W_{0}$ and type of mass fallback is depicted. As it was already pointed out by Weinberg (1993), models on the left of the vertical black line (so called ’Weinberg’s cliff’) are dissolving on the dynamical time-scale because of the strong mass loss connected with stellar evolution. Models on the right of this line evolve on the relaxation time scale, and usually before dissolution they undergo core collapse. For the Mocca models this is true with the exception of the tidally filling model with $W_{0}=6$ and with mass fallback ON. For this model, despite the evolution time-scale being proportional to the relaxation time and the relatively early core collapse, it enters an abrupt dissolution phase and very quickly dissolves, still being on the right side of the ’Weinberg’s cliff’. As we argued above, this dissolution is connected with the formation and evolution of a BHS in strongly tidally striped clusters. This model is an example of the third dissolution mechanism.

3.2 Cluster dissolution - black hole subsystem properties and evolution

In the previous Section we showed that the BHS seems to play a crucial role in the dissolution of the tidally filling star clusters. Now, we will try to understand why that is by looking at the BHS evolution and its properties. It is well known from the theoretical paper by Breen & Heggie (2013) and simulation results discussed by Arca Sedda et al. (2018), Askar et al. (2018) and Kremer et al. (2019) that the evolution of the BHSs depends strongly on the cluster half-mass relaxation time and it is connected with the continuous removal of BHs in dynamical interactions - first most massive and then less massive ones. Such a decreasing number of BHs is accompanied with smaller BHS size, larger BHS characteristic density, smaller average BH masses, higher fraction of BHs in binaries and larger average mass of ordinary stars mixed with BHs in the BHS. Using these relations we will try to understand the reason for tidally filling cluster dissolution by BHSs.

Figures 6 show the evolution of the overall Lagrangian radii and BHS Lagrangian radii for tidally filling and tidally under-filling models. From the inspection of the Figure we can clearly see the initial overall cluster expansion powered by the stellar evolution mass loss. For the tidally under-filling model, expansion is stronger than in the tidally filling model, because for such a model, the cluster can freely expand until it reaches the tidal radius. Then, for both models we can observe the mass segregation of BHs. For the tidally filling model it stops for the innermost $10\%$ Lagrangian radius at about 0.5 for the fraction of cluster bound mass, and for the tidally under-filling model at about 0.75 for the fraction of cluster bound mass. The mass segregation for the other BH Lagrangian radii continues, and it stops at about 0.3, or never stops for tidally filling and tidally under-filling models, respectively. For the tidally filling model we can clearly see that BH mass segregation is connected with the expansion of the central parts of the cluster (except $1\%$ overall Lagrangian radius, which is dominated by BHs) and decreasing outer overall Lagrangian radii (see the $50\%$ Lagrangian radius), which is connected with strong tidal striping of the cluster. This behaviour, as we will see later, seems to be very important for understanding the final abrupt cluster dissolution. For the tidally under-filling model, the inner overall Lagrangian radii starts to decrease when the BHS contracts to support cluster structure and needs to generate more energy (decrease of the $1\%$ BHS Lagrangian radius). Half-mass radius is still increasing in accord with central energy generation and the tidal stripping is not too strong.

The described behaviour above can be easily understood on the basis of the theoretical arguments provided by Breen & Heggie (2013) and results of simulations discussed by Arca Sedda et al. (2018); Askar et al. (2018). The first part of the balanced evolution of the BHS is characterized by a nearly constant number of BHs and a nearly constant ratio between BHS mass and the cluster mass. The duration of this phase depends on the time-scale of BH mass segregation. When the most massive BHs are kicked out in dynamical interactions the mass of the BHS is decreasing together with the average BH mass (see Figures 7 and 8). The duration of this phase of evolution depends on the half-mass relaxation time. The shorter the half-mass relaxation time, the faster the BHS evolution, the faster the decrease of the average BH mass, and the faster the increase of the BHS density (smaller innermost Lagrangian radius). This general picture is slightly modified for the tidally filling model for the fraction of cluster bound mass smaller than 0.4 and larger than 0.2, when the mass ratio between BHS mass and the total cluster mass is nearly constant instead of decreasing. Such departure from the standard theoretical picture (Breen & Heggie, 2013) can be attributed to the strongly decreasing total cluster mass due to tidal striping. The required energy flow through the half-mass radius is constantly decreasing, so the BHS evolves slower and its mass is changing nearly in accord with the cluster mass. There is some decrease of the average BH mass, due to the removal of the most massive BHs in dynamical interactions. At the point when the the fraction of cluster bound mass is equal to about 0.2, the BHS starts to decouple from the rest of the system. The average BH mass is nearly constant, which suggests a very small rate of energy generation by BHS in dynamical interactions. The BHS Lagrangian radii starts to slowly increase and the cluster stars are quickly tidally stripped, which suggests a loss of dynamical equilibrium (see in Figure 7 the steep increase of the ratio between BHS mass and the cluster mass).

Already Fukushige & Heggie (1995); Whitehead et al. (2013) and Contenta et al. (2015) showed that the abrupt cluster dissolution powered by strong mass loss due to stellar evolution is connected with the overall loss of cluster dynamical equilibrium. Banerjee & Kroupa (2011) in turn showed, for small N clusters, that so the called ’dark clusters’ are formed from the BHS, which forces the rest of the system to abrupt dissolution, because of the loss of dynamical equilibrium for the luminous objects (see their Figure 1, top panel). Figure 9 shows the virial ratios computed for the whole system, for BHS and other objects, for the tidally filling model (top panel) and the tidally under-filling model (bottom panel), respectively. In order to compute virial ratios for BHS and other objects from the Mocca models we needed to use the detailed snapshots which are less frequently stored (every 100 Myr) than the standard output (every few Myr). Therefore, there are larger gaps between points for virial ratios computed from the snapshots. From Figure 9 (top panel) we can clearly see that the cluster dissolution for tidally filling model is indeed connected with the decoupling of the BHS from the rest of the cluster and the loss of dynamical equilibrium by other objects. Unfortunately, because of the scarcity of snapshots very close to the cluster dissolution time, we cannot observe the final values of the virial ratios for all cluster components. We also clearly see, that the BHS stops to collapse further at around the fraction of cluster bound mass equal to about 0.2, which coincidences with the slow increase of the BHS Lagrangian radii. The BHS starts to disrupt itself. Its evolution follows the well known behaviour of low-N systems. For the tidally under-filling cluster (bottom panel) we clearly see mass segregation of the BHs and formation of a BHS, but at the fraction of cluster bound mass equal to about 0.65, further decrease of the BHS virial ratio is stopped. It seems that this is connected with the faster removal of the most massive BHs in dynamical interactions than for tidally filling models, and is accompanied by the decrease of BHS Lagrangian radii. The large fluctuations observed for the BHS virial ratios are connected with the decreasing number of BHs. Figure 9 clearly confirms findings by Banerjee & Kroupa (2011) and shows that a ’dark cluster’ can be formed also for star clusters of the size of globular clusters, not only for open clusters.

3.3 Cluster dissolution - Global similarities

Breen & Heggie (2013) in their model approximated a star cluster harbouring a BHS by a simple two-component model. They assumed, following the Hénon principle (Hénon, 1975), that the rate of energy generation by BHS is regulated by the energy demand of the bulk of the system, which is usually thought of as the energy flux at the half-mass radius. Using this approximation, they derived a formula, which relates the half-mass radius of the BHS to the overall cluster half-mass radius. Using equation 4 from Breen & Heggie (2013) and the parameters of the BHS and the clusters from the Mocca simulations we plotted on Figure 10 the ratio between the half-mass radius of the BHS to the overall cluster half-mass radius for Mocca models, characterized by a different: number of objects, binary fractions, tidal radii and $W_{0}$. The model parameters are described in the Figure caption. Except models with $W_{0}=3$ and $W_{0}=9$, all other models are grouped in two separate families. First, connected with all models having initial binary fraction equal to $95\%$, regardless of the number of objects or $R_{t}$. Second, connected with models having small binary fractions, also regardless of the number of objects and $R_{t}$. As it was already discussed before, the model with $W_{0}=3$ evolves and dissolves very quickly, so the BHs are unable to mass segregate, and therefore the ratio between the BHS half-mass radius and overall half-mas radius is nearly constant. Model with $W_{0}=9$ is initially very dense and has very short half-mass relaxation time, so the BHS is formed and burned out quickly, therefore the ratio between half-mass radii is decreasing. The observed differences between models with a small and large binary fraction are connected with the fact that small binary fraction models are initially able to retain a larger fraction of BHs on the cluster unit mass. This rather unexpected conclusion (taking into account the fact that models with a very large binary fraction contain nearly twice the number of stars compared to models with a low binary fraction) is connected with the fact of very different initial binary properties in both models. Models with a large binary fraction follow period and mass ratio distributions according to Kroupa (1995). For such models in the Mocca-Survey Database, the mass ratio is practically one for MS stars which can form BHs, and binary periods can be very large. For such binary parameters most of the wide binaries will be disrupted, even for the massive MS stars, by two successive (practically at the same time) SNe natal kicks. A significant fraction of such newly formed BHs will escape from the system, except the most massive ones. This process substantially decreases the number of BHs in binaries and increases the number of single BHs kept in the system. Also, binaries that contain two massive stars and have small initial separations can merge during the common envelope evolution, producing a single BH instead of a BBH. Therefore, the number of single BHs will be much larger than the number following from the initial number of single stars, but much smaller than the number following from the initial total number of stars (single stars and binary components). Contrarily, for models with initially low binary fractions, there are fewer binaries to begin with and they are rather hard and will usually merge before the first SNe and form one massive BH with a small natal kick. Therefore, they will be kept in the system substantially increasing the number of formed (and kept) BHs from single stars. Accordingly, the number of BHs in binaries is very low - most BHs are single. According to equation 4 in Breen & Heggie (2013) the larger the BH mass fraction the larger the ratio of the half-mass radii for BHs and the whole system.

Figure 11 shows the evolution of the ratios between half-mass radii of BHS and the whole system, similar to Figure 10, but this time we are showing the half-mass radii from the Mocca models and not the evaluated values from equation 4 in Breen & Heggie (2013). Generally, both figures are very similar, after the initial phases of strong stellar evolution mass loss and BH mass segregation, which leads to the balanced cluster evolution. Models with lower numbers of objects achieve the phase of the balanced evolution for larger values of the fraction of cluster bound mass. This behaviour is connected with much faster BH mass segregation for smaller systems than for larger ones. The models with smaller binary fraction also show larger half-mass radii ratios as on Figure 10, but this time the difference is less pronounced. This suggests, as expected, that the structure of the BHSs in star clusters is more complicated than a simple two-component model. It also suggests, that the simple Hénon principle (Hénon, 1975) properly captures the global evolution of star clusters, which have a complicated structure and contain a mixture of different kinds of objects responsible for the central energy generation. It is worth pointing out that the levels of the ratios of BHS half-mass radius and cluster half-mass radius are very similar on Figures 10 and 11, which further strongly supports the theoretical picture of BHSs given by Breen & Heggie (2013).

The evolution of the ratio between $10\%$ BHS Lagrangian radius and cluster half-mass radius for different models is shown on Figure 12. Generally, the main features are very similar to Figure 11, but now the differences between models with low and large binary fractions are practically invisible. All models during balanced evolution have very similar ratios between $10\%$ BHS Lagrangian radius and cluster half-mass radius. This suggests that the energy generation by BHS is confined in the inner $10\%$ BHS Lagrangian radius and provides a universal path for cluster evolution regardless of the global cluster parameters.

3.4 Cluster dissolution - Summary of the most important features

In the previous Sections, we have provided arguments that tidally filling clusters harbouring BHSs in the centres can abruptly dissolve, provided that the clusters are described initially by a moderate King concentration parameter (around $W_{0}=6$), and cluster masses before the dissolution are about $10\%$ of the initial cluster masses. In this Section we will provide a short and coherent picture to take away from the paper.

The evolution of a tidally filling cluster is controlled by the following physical processes (if not stated otherwise all points below describe evolution of tidally filling models):

initial mass loss is connected with stellar/binary evolution. All models regardless of their initial properties and for the canonical IMF (Kroupa, 2001) are controlled by the mass loss up to the time of around 200 Myr and the fraction of cluster bound mass equal to about 0.7 - see Figures 1 and 2; 2. 2.

then relaxation process takes over and BHs start to segregate and form a BHS. It seems that the BHS mass segregation stops around mass fraction equal to 0.45, depending on the number of stars, and takes slightly longer for the outer BHS Lagrangian radii. For tidally under-filling clusters, BH mass segregation is much faster and never stops for the outer BHS Lagrangian radii. This phase of evolution is characterized by: slightly decreasing average BH mass, constant BHS mass to cluster mass ratio and nearly constant Lagrangian radii for the BHS - see Figures 8, 7 and 6 top panel, respectively. Massive BBHs generate enough energy to support cluster evolution without strongly changing BHS properties; 3. 3.

a bit later, the fraction of cluster bound mass for a tidally filling cluster starts to differ from that for a tidally under-filling cluster - see Figure 1. This happens at mass fraction of about 0.3. The tidally under-filling cluster at this moment does not contain any BHS and energy is generated by stellar-type binaries. During this phase the tidal stripping becomes more and more efficient and the overall half-mass radius starts to slowly decrease, but BHS Lagrangian radii are still nearly constant. Accordingly, BHS mass slowly increases compared to the cluster mass and the average BH mass is slightly decreasing - see Figures 6, 7 and 8, respectively; 4. 4.

at around mass fraction equal to 0.2, the cluster enters the phase of abrupt dissolution. During this phase, the average BH mass is nearly constant and the mass ratio between BHS mass and cluster mass starts to strongly increase. BHS Lagrangian radii start to slowly increase and the overall Lagrangian radii start to decrease with larger and larger speed. The virial ratio for BHS stops to decrease and becomes nearly constant, but the virial ratios for objects other than BHs start to increase showing enhanced removal of stars. The Cluster starts to evolve on a dynamical time-scale, not on a relaxation time-scale (Banerjee & Kroupa, 2011) - see Figures 8, 7, 6 top panel and 9, respectively; 5. 5.

the strong energy generation by massive BHS, and as a consequence enhanced tidal stripping, are the main mechanisms responsible for abrupt cluster dissolution. The reasons are following:

1. (a)

as it was pointed out by Arca Sedda et al. (2018) and Askar et al. (2018) BHs in BHS are very well mixed with other stars. The energy transport from BHS to the cluster halo occurs via stars kicked out from the BHS. Massive BBHs are extremely hard from the point of view of other, less massive, objects in the BHS, so during dynamical interactions, stars, and from time to time also BHs, are vigorously kicked out from the BHS; 2. (b)

some objects kicked out from the BHS will also escape from the cluster. The number of escaped objects will depend on the central escape velocity; 3. (c)

cluster is constantly tidally stripping, but the internal cluster structure is not strongly changing (see Figures 6 top panel and 8), so the escape velocity becomes smaller and smaller (see Figure 13 top panel) and more and more objects can be kicked out from the cluster - see Figure 13 bottom panel; 4. (d)

the positive feedback between strong energy generation, high flow of escaping stars, decreasing escape velocity and enhanced tidal stripping leads finally to loss of cluster dynamical equilibrium and abrupt dissolution - see Figure 13 bottom panel.

In our opinion, the above described third mechanism of cluster dissolution is universal and should work for all clusters, which harbour BHSs to the end of their life, and are exposed to a strong tidal field, regardless of the initial cluster structure or surrounding environment.

4 Discussions and Conclusions

In the previous Section we showed that tidally filling clusters with BHS in the center can show abrupt dissolution, provided that the cluster is described initially by a moderate King concentration parameter (around $W_{0}=6$), and the cluster mass before dissolution is about $10\%$ of the initial cluster mass. The dissolution is controlled by a strong energy generation in dynamical interactions inside the BHS and strong tidal stripping which lead to sudden loss of dynamical equilibrium. When the cluster loses a significant fraction of its initial mass, the BHS decouples from the rest of the cluster. The outer parts of the system cannot accommodate energy flux, because of increasing tidal stripping. Models that have larger $W_{0}$ or are strongly tidally under-filling are not able to keep their BHSs until the late phases of cluster evolution. They are characterized by slow dissolution controlled by the relaxation process and usually are core collapsed clusters. Tidally filling models with $W_{0}$ smaller than about $3$ will show abrupt dissolution which is caused by the strong early mass loss connected with stellar/binary evolution and extremely fast tidal stripping. They usually do not enter the core collapse phase.

The above picture of the abrupt cluster dissolution by BHS energy generation and strong tidal field influence suggests that such a mechanism should operate not only in tidally filling clusters. For tidally under-filling clusters and/or clusters with larger $W_{0}$ we need initially a physical process which will quickly and strongly expand the cluster and/or decrease high initial concentration. Actually, such a process is operating in the initial phases of the cluster evolution, namely mass loss connected with the stellar/binary evolution. Only we need to make it more efficient. This can be achieved by making the IMF top-heavy. To check this hypothesis we started a new project to investigate how the cluster dissolution process depends on the IMF power-law index for massive stars (Wang et al., 2019). The very preliminary results are shown on Figure 14 for a model with total mass equal to $500000$ $M_{\odot}$, $W_{0}=6$, strongly tidally under-filling with $R_{t}/R_{h}=80$, binary fraction equal to $90\%$, fallback ON and different IMF power-law index $\alpha$. As expected, models with top-heavy IMF show an abrupt dissolution which is very similar to the one described in this paper for tidally filling clusters with BHS. Models with extremely low $\alpha$ show dissolution which is powered by the initial mass loss connected with rapid stellar evolution not with energy generation by BHs. Very low $\alpha$ means a very large number of very massive stars which will very quickly lose most of their mass due to stellar evolution. Some indications about abrupt dissolution by top-heavy IMF were also provided by Chatterjee et al. (2017), but the authors did not discuss in detail the physical processes leading to such cluster dissolution. They only mentioned that strong energy generation could be responsible for such behaviour.

There are more and more observational evidences (e.g. Haghi et al., 2017; Jeřábková et al., 2018; Kalari et al., 2018; Schneider et al., 2018; Hosek et al., 2019) and theory arguments (e.g. Marks et al., 2012) that star clusters are born with top-heavy IMF. As it was pointed out by Marks et al. (2012), how strongly IMF is top-heavy depends on the environment properties (density, mass) of stellar forming regions. If globular clusters were formed in an environment with strong gas ram pressure, a stripping and quickly changing galactic tidal field (e.g. Kruijssen, 2014; Kruijssen et al., 2018; Renaud, 2018) we should expect that globular clusters that survive up to the Hubble time should be massive with large central concentrations. Probably clusters with smaller galactocentric distances were tidally filling or not strongly tidally under-filling. Cluster with larger distances, probably migrated out and should be strongly tidally under-filling (Kruijssen, 2014). Further spatial evolution of globular clusters will be shaped. among others, by two physical processes. One connected with dynamical friction, which drag clusters to the galactic center, and other connected with cluster dissolution, in some of them powered by BHSs. In both processes, globular clusters will be lost relatively early from the observational point of view. The latter process will lead to the buildup of the nuclear star cluster, removal of the most massive clusters and is the more effective one from the point of view of the number of retained BHs and the possible BH-BH mergers. The former process will relatively quickly disrupt clusters leading to a strong increase of BHs and BBHs populations in the galactic halo/bulge. Such a large population of BHs should be accompanied with some observational features connected with microlensing BHs (e.g. Rybicki et al., 2018) and provide more firm constrains on the number of BHs floating in the Galactic halo.

According to the third mechanism of cluster dissolution the cluster dispersion is abrupt in time. The final stages of such a cluster should contain in the core a dark component surrounded by luminous stars. The dark core should measure about 1 pc and the cluster half-mass radius should be around 10 pc. The possibility of the formation of such ’dark clusters’ in the population of open clusters has been predicted by Banerjee & Kroupa (2011), but in this paper we extend this to the population of globular clusters. From the observational point of view, such clusters should be rather difficult to observe because they are very short living, contrary to the other type of ’dark clusters’ proposed by Askar et al. (2017b). Such ’dark clusters’ are composed of massive intermediate mass BH surrounded by bright stars with a smaller total mass. Such clusters are characterized by much larger central velocity dispersion and can live in such a stage much longer than clusters described in this paper. It seems that the two types of ’dark clusters’ proposed in the literature have very different observational properties and should not be confused observationally.

The Mocca code used to create the Mocca Survey Database I does not include post-Newtonian effects, which could be important during chaotic, strong dynamical interactions between BBHs and BHs in a dense and large BHSs. Because evolution of such systems is discussed in this paper, therefore we decided to comment on the possible influence of such effects on the cluster global structure and evolution. There are two main effects connected with gravitational wave emission: kinetic energy dissipation leading to an enhanced merging rate, and merger velocity kicks leading to an enhanced escape rate. Both effects will lead to the decrease of the BH number in the system. According to discussions in the previous Sections, it was argued that the most massive BBHs, composed of nearly equal mass BHs, are responsible for most of the energy generated in BHS and for slowing down the cluster evolution rate. If such a BBH will merge (because of gravitational wave radiation) they will be probably kicked out from the system due to higher recoil kicks for high mass ratio BBHs (see e.g. Morawski et al., 2018). Earlier removal of such a binary will lead to the formation of another, a bit less massive BBH, which in order to provide energy needed to support cluster structure will have to more frequently dynamically interact with other objects. This process, in turn will lead to the increase of cluster central density, faster dissolution of BHS and faster cluster evolution, preventing, in a most extreme situation, abrupt cluster dissolution. Taking into account strong stochasticity of dynamical interactions, we believe that such a scenario is rather unlikely. Firstly, most of the energy provided by massive BBH is generated with interactions with low mass ordinary stars, which are responsible for energy transport in the cluster. Such interactions will slowly increase the hardness of BBH, finally leading to its escape. Secondly, mergers between interacting BHs require very close approaches (comparable with the BH Schwarzschild radius), which in turn mean very hard BBH. It is possible that such hard BBH will be kicked out from the system before it can be prone for gravitational wave radiation mergers. Summarizing, we believe that the scenario of abrupt cluster dissolution discussed in this paper will be only slightly, if at all, affected by adding the post-Newtonian effects.

In this paper, we have presented on the basis of Mocca simulations of real star cluster evolution, strong evidence of the existence of a third mechanism leading to cluster dissolution. In addition to cluster dissolution connected with strong mass loss due to stellar/binary evolution or connected with mass loss controlled by the relaxation process, we propose a dissolution mechanism connected with the presence of BHSs in tidally filling clusters. Energy generation in BHS together with strong tidal stripping leads to abrupt cluster dissolution. The dissolution process is connected with the loss of cluster dynamical equilibrium, similarly as it is in the case of dissolution powered by stellar evolution mass loss. We have argued that the third mechanism should be universal and will occur for star clusters embedded in a strong tidal field, which harbour massive BHS until the last phases of their evolution. Therefore, it should also work for a tidally under-filling cluster with top-heavy IMFs. This possibility will be further investigated in detail in the next paper.

Acknowledgements

We wish to thank Douglas Heggie for the extremely insightful comments and in-depth discussions. MG and AL were partially supported by the Polish National Science Center (NCN) through the grant UMO-2016/23/B/ST9/02732. AA is currently supported by the Carl Tryggers Foundation for Scientific Research through the grant CTS 17:113 and was partially supported by NCN, Poland, through the grants UMO-2016/23/B/ST9/02732. LW in this work was supported by the funding from Alexander von Humboldt Foundation. AH was supported by Polish National Science Center grant 2016/20/S/ST9/00162. RS has been supported by the Alexander-von-Humboldt Polish Honorary Research Fellowship of the Foundation for Polish Science, and by the National Natural Science Foundation of China (NSFC), grant 11673032. We acknowledge support from the Chinese Academy of Sciences through the Silk Road Project at the National Astronomical Observatories, Chinese Academy of Sciences (NAOC), through the ’Qianren’ special foreign experts programme of China, and through the Sino-German collaboration program GZ1284. This work benefited from support by the International Space Science Institute (ISSI), Bern, Switzerland, through its International Team programme ref. no. 393 The Evolution of Rich Stellar Populations & BH Binaries (2017-18).