# The survival of star clusters with black hole subsystems

Long Wang1,2,3
1 Argelander Institut Für Astronomie, Auf Dem Hügel 71, 53121, Bonn, Germany
2 Helmholtz-Institut für Strahlen- und Kernphysik, University of Bonn,
Nussallee 14-16, D-53115 Bonn, Germany
3 RIKEN Center for Computational Science, 7-1-26 Minatojima-minami-machi,
Chuo-ku, Kobe, Hyogo 650-0047, Japan
E-mail<EMAIL_ADDRESS>
(Accepted –. Received –; in original form –)

###### Abstract

Recent observations have detected top-heavy IMFs in dense star forming regions
like the Arches cluster. Whether such IMFs also exist in old dense stellar
systems like globular clusters is difficult to constrain, because massive
stars already became black holes (BHs) and neutron stars (NSs). However,
studies of stellar dynamics find that BHs/NSs influence the long-term
evolution of star clusters. Following Breen & Heggie (2013) and by carrying
out two-component $N$-body simulations, we demonstrate how this dynamical
impact connects with the shape of IMFs. By investigating the energy balance
between the BH subsystem and the global, we find that to properly describe the
evolution of clusters, a corrected two-body relaxation time,
$T_{\mathrm{rh,p}}=T_{\mathrm{rh}}/\psi$, is necessary. Because $\psi$ depends
on the total mass fraction of BHs, $M_{\mathrm{2}}/M$, and the mass ratio,
$m_{\mathrm{2}}/m_{\mathrm{1}}$, the cluster dissolution time is sensitive to
the property of BHs or IMFs. Especially, the escape rate of BHs via ejections
from few-body encounters are linked to mass segregation. In strong tidal
fields, top-heavy IMFs easily lead to the fast dissolution of star clusters
and the formation of BH-dominant dark clusters, which suggests that the
observed massive GCs with dense cores are unlikely to have extreme top-heavy
IMFs. With the future observations of gravitational waves providing unique
information of BHs/NSs, it is possible to combine the multi-message
observations to have better constrains on the IMFs of old star clusters.

###### keywords:

methods: numerical – galaxies: star clusters: general – stars: black holes.

††pagerange: The survival of star clusters with black hole subsystems–The
survival of star clusters with black hole subsystems††pubyear: 2019

##  1 Introduction

Whether the initial mass function (IMF) is environmentally dependent is a long
existing question. The observations of present-day star forming regions (young
star clusters and associations) in the Galaxy and Magellanic clouds show no
strong evidence that the high-mass end of IMFs ($>1M_{\odot}$) depends on the
density for the range of $\rho_{\mathrm{0}}=35-3\times
10^{4}M_{\odot}\mathrm{pc}^{-3}$, where $\rho_{\mathrm{0}}$ is the central
density of clusters, and the metallicity $Z=0.2-1Z_{\odot}$ (Bastian, Covey &
Meyer, 2010). However, extremely dense environments like the Arches cluster
($\rho\approx 2\times 10^{5}M_{\odot}\mathrm{pc}^{-3}$), the Galactic center
and 30 Doradus are found to contain top-heavy IMFs with heavy-end exponent
$\alpha\approx-1.7\sim-1.9$ (Lu et al., 2013; Schneider, et al., 2018; Hosek,
et al., 2019). Recently, by measuring the 13CO/C18O ratio in starburst
galaxies at high redshift ($z\sim 2-3$), Zhang et al. (2018) found that a top-
heavy IMF with $\alpha\approx-2.1$ is necessary to explain the star formation
rate there. Thus, there is a possible trend that IMF might be top-heavy in
high-density environments. If this is the case, an interesting question is
whether old and dense globular clusters (GCs) in the Galaxy also contain top-
heavy IMFs. It is difficult to determine IMFs in GCs since all massive stars
have already evolved to black holes (BHs) and neutron stars (NSs) a long time
ago. A semi-analytic work by Marks et al. (2012) estimated the IMFs of GCs and
suggested that GCs can contain extreme top-heavy IMFs
($\alpha\approx-0.7\sim-2.3$) depending on density and metallicity. However, a
recent study by Baumgardt & Sollima (2017) compared the scaled $N$-body models
with observational data of $35$ GCs found that the present-day mass function
of GCs agrees well with Kroupa (2001); Chabrier (2003) IMFs. Thus, the studies
of IMFs in GCs are controversial and there is no star-by-star $N$-body
simulations of GCs yet to clarify this due to the time-consuming computing.

One significant influence of IMFs (heavy-mass end) is the number fraction of
BHs and NSs formed in star clusters. By including the BH subsystems in fast
Monte-Carlo simulations, Chatterjee, Rodriguez, & Rasio (2017) and Giersz, et
al. (2019) found that tidal filling GCs with top-heavy IMFs dissolve much
faster than the models with Kroupa/Chabrier IMFs, thus more difficult to
survive until today. Therefore, the dynamical effect of BHs may provide a
strong constraint on the IMFs of star clusters, especially in GCs.

It is very useful to develop a theory how BHs dynamically affect the evolution
of star clusters, in order to provide a fast method to describe the properties
of BHs and star clusters without expensive $N$-body models. Hénon (1961, 1975)
showed that the long-term evolution of star clusters are guided by the balance
of energy flux, i.e. the energy generation in the cluster core is regulated by
the energy demand of the system. Based on this idea, Breen & Heggie (2013)
suggested that in star clusters with many BHs, the energy generation is
dominated by the few-body encounters between BH binaries and singles in a
centrally concentrated BH subsystem. Thus, the balance of energy flux can be
described by (Eq. 1 in Breen & Heggie, 2013)

$\frac{|E|}{T_{\mathrm{rh}}}\approx
k\frac{|E_{\mathrm{2}}|}{T_{\mathrm{rh,2}}},$ (1)

where $|E_{\mathrm{2}}|$ and $|E|$ are the total energy of the BH subsystem
and the global system respectively; $T_{\mathrm{rh}}$ and $T_{\mathrm{rh,2}}$
are the two-body relaxation times measured at the half-mass radius of the
system, $R_{\mathrm{h}}$, and the BH subsystem, $R_{\mathrm{h,2}}$,
respectively. Hereafter, the suffix “2” and “1” represent BH subsystems and
light star components separately. When the global parameters of a cluster are
known, this relation provides a constrain on the properties of BH subsystems.
It results in several important features of the long-term evolution of star
clusters:

  * •

During the balanced stage, $R_{\mathrm{h}}$ increases due to the heating of
BHs and $R_{\mathrm{h,2}}/R_{\mathrm{h}}$ has a certain relation (Eq. 4 in
Breen & Heggie, 2013).

  * •

BHs are trapped in the deep potential of the cluster, thus the escape of BHs
via tidal stripping can be ignored. The ejection of BHs via few-body
encounters between BH binaries and singles in the cluster core is the major
channel to remove BHs from the cluster. Since the energy generation rate via
the few-body encounters is controlled by the global energy requirement, the
escape rate of BHs is independent of their relaxation time,
$T_{\mathrm{rh,2}}$, but scales with $T_{\mathrm{rh}}$ of the global system
(Eq. 9 in Breen & Heggie, 2013).

  * •

In tidally limited systems, if the total mass fraction of BHs is small, i.e.
$M_{\mathrm{2}}/M<0.1$, the BHs are depleted faster and $M_{\mathrm{2}}/M$
decreases as clusters evolve. Otherwise $M_{\mathrm{2}}/M$ increases and
finally “dark clusters” form (Banerjee & Kroupa, 2011; Giersz, et al., 2019).

The first two points assume that BHs are not the dominant component in star
clusters, which is the general case for canonical IMFs (Kroupa, 2001;
Chabrier, 2003).

On the other hand, with the help of $N$-body simulations, Baumgardt (2001);
Baumgardt & Makino (2003) developed a theory for the mass loss of star
clusters and showed that the dissolution time of the systems is (by putting
Eq. 1 to Eq. 5 from Baumgardt & Makino, 2003, and canceling
$R_{\mathrm{G}}/V_{\mathrm{G}}$ in Eq. 5)

$T_{\mathrm{dis}}\sim
T_{\mathrm{rh}}^{x}T_{\mathrm{cr}}^{1-x}(R_{\mathrm{t}}/R_{\mathrm{h}})^{3/2},$
(2)

where $T_{\mathrm{cr}}$ is the crossing time, $R_{\mathrm{t}}$ is the tidal
radius of the system. The power index $x$ depends on the structure of the
cluster and the tidal field, and is around $0.75$. The term
$T_{\mathrm{cr}}^{1-x}$ represents the back-scatter effect of potential
escapers and the term $(R_{\mathrm{t}}/R_{\mathrm{h}})^{3/2}$ is a radius
scaling factor that scales the half-mass radius by the tidal radius. This
theory provides a simple solution to evolve the masses of star clusters by
given initial conditions. However, the effect of BH subsystems is not included
in their study, thus it is unknown whether the theory can be extended to
clusters with top-heavy IMFs.

In this work, as a follow-up study of Breen & Heggie (2013), we extend the
energy-balance theory to the region of high $M_{\mathrm{2}}/M$, which
represents star clusters with top-heavy IMFs. In Section 2, how top-heavy IMFs
affect $M_{\mathrm{2}}/M$ is shown. Then we describe the initial conditions of
the two-component $N$-body models in Section 3. We provide the results and the
analysis about the energy-flux balance, the mass loss of BHs and stars in
isolated environments in Section 4 and in the tidally filling case in Section
5. Finally, we discuss our results and draw conclusions in Section 6.

##  2 Top-heavy IMF

The top-heavy IMFs can be defined by varying the slope, $\alpha_{\mathrm{3}}$,
in the heavy end of IMFs ($1-150M_{\odot}$) based on the multi-component
power-law IMF (Kroupa, 2001),

$\displaystyle\xi(m)\propto m^{\alpha_{\mathrm{i}}}$ (3)
$\displaystyle\alpha_{\mathrm{1}}$ $\displaystyle=-1.3,$ $\displaystyle
0.08\leq m/M_{\odot}$ $\displaystyle<0.50$ $\displaystyle\alpha_{\mathrm{2}}$
$\displaystyle=-2.3,$ $\displaystyle 0.50\leq m/M_{\odot}$
$\displaystyle<1.00$ $\displaystyle\alpha_{\mathrm{3}}$ , $\displaystyle
1.00\leq m/M_{\odot}$ $\displaystyle<150.0.$

The total mass fraction of BHs, $M_{\mathrm{2}}/M$, after $100$ Myr evolution
can be estimated by using the updated sse/bse code (Hurley, Pols, & Tout,
2000; Banerjee, et al., 2019). The models of stellar winds and formation of
NSs and BHs are detailed discussed in Banerjee, et al. (2019). The remnant
mass is based on the rapid supernovae scenario from Fryer, et al. (2012) and
the calculation also includes the pair-instability supernovae (Belczynski, et
al., 2016). The result is shown in Table 1. Here all formed BHs are counted in
the calculation and the metallicity of stars is $Z=0.001$. A part of BHs may
immediately escape the clusters due to high kick velocities after the
supernovae. This effect is ignored in our analysis. For the case of a
canonical IMF of Kroupa (2001) with $\alpha_{\mathrm{3}}=-2.3$, after $100$
Myr, BHs contribute about $7\%$ of the total mass and about $20\%$ of the
initial mass is lost via stellar evolution (winds) and supernovae. When
$\alpha_{\mathrm{3}}=-1.7$, which is the case for the Arches cluster (Hosek,
et al., 2019) and the Galactic center (Lu et al., 2013), the mass fraction of
BHs grows significantly (to $38\%$). The average mass ratio for BHs and stars,
$\langle m_{\mathrm{2}}\rangle/\langle m_{\mathrm{1}}\rangle$, is about $40$
for metallicity $Z=0.001$. This result can be altered (not significantly) by
choosing different assumptions in the stellar evolution model. But in general
the shape of the IMF changes $M_{\mathrm{2}}/M$ strongly.

Table 1: The BH properties after $100$ Myr stellar evolution for different IMFs with the metallicity, $Z=0.001$. Columns show the mass fraction of BHs, $M_{\mathrm{2}}/M$, the stellar wind mass loss normalized to the initial total mass, $M_{se}/M(0)$, the number fraction of BHs, $N_{\mathrm{2}}/N$, and the average mass ratio for BHs and stars, $\langle m_{\mathrm{2}}\rangle/\langle m_{\mathrm{1}}\rangle$. $\alpha_{\mathrm{3}}$ | $M_{\mathrm{2}}/M$ | $M_{se}/M(0)$ | $N_{\mathrm{2}}/N$ | $\langle m_{\mathrm{2}}\rangle/\langle m_{\mathrm{1}}\rangle$
---|---|---|---|---
-1.5 | 0.554 | 0.536 | 0.0319 | 37.8
-1.7 | 0.376 | 0.466 | 0.0151 | 39.3
-2.0 | 0.182 | 0.334 | 0.00522 | 42.5
-2.3 | 0.0733 | 0.204 | 0.00181 | 43.7

##  3 N-body models

A BH is much more massive than a light star. It is expected that the stellar
dynamical effect of BHs play an important role for the long-term evolution of
star clusters. With top-heavy IMFs, this effect is more pronounced. In order
to investigate this, we carry out a series of $N$-body simulations of star
clusters with two components and no stellar evolution (similar to Breen &
Heggie, 2013). The spherically symmetric Plummer models are employed for
initial conditions. Two components are homogenously mixed initially with the
same density and velocity profile. Each component has an equal mass for all
members. The heavy and light components represent BHs and light stars
separately. Several physical processes involved in the evolution of a real
star cluster are ignored in these simple models. However, in the first step,
it is important to obtain a better theoretical understanding by investigating
the simplest model with the interesting physical process well isolated.

The initial number fraction of the heavy component, $N_{\mathrm{2}}/N(0)$, and
the average stellar mass ratio, $m_{\mathrm{2}}/m_{\mathrm{1}}$, are varied
for different models. For convenience, hereafter “m[value]” and “N[value]”
represent different $m_{\mathrm{2}}/m_{\mathrm{1}}$ and $N_{\mathrm{2}}/N(0)$
separately (e.g. m2.5 is $m_{\mathrm{2}}/m_{\mathrm{1}}=2.5$ and N0.01 is
$N_{\mathrm{2}}/N(0)=0.01$). The initial total mass ratio of two components,
$M_{\mathrm{2}}/M(0)$, are listed in Table 2. $M_{\mathrm{2}}/M(0)$ vary from
$2.5\%$ to $67\%$, thus this model set represents a wide range of top-heavy
IMF. $m_{\mathrm{2}}/m_{\mathrm{1}}$ includes the values of $2.5$ to $40$,
corresponding to the masses of NSs and BHs. The initial total number of
objects have two values, $32k$ and $64k$, as shown in Table 2.

For each model, we consider the star clusters in an isolated environment and
in a point-mass potential with a circular orbit. A star is defined as an
escaper if its energy is positive and its distance to the cluster center is
larger than $R_{\mathrm{esc}}$. In the isolated case,
$R_{\mathrm{esc}}=20~{}R_{\mathrm{h}}$. In the tidal filling case,
$R_{\mathrm{esc}}=2~{}R_{\mathrm{tid}}$, and initially,
$R_{\mathrm{tid}}/R_{\mathrm{h}}\approx 7.3$, for all models, which is close
to the case of the tidal filling King model with $W=6\sim 7$.

The nbody6++gpu code (Wang et al., 2015) is used to perform the $N$-body
simulations. The code is designed to simulate star clusters with
regularization algorithms to accurately deal with the strong interaction
between stars and BHs, and therefore fits the purpose of this study (Aarseth,
2003, and reference therein).

Table 2: The initial $M_{\mathrm{2}}/M(0)$ for all two-component $N$-body models. Rows show the same $m_{\mathrm{2}}/m_{\mathrm{1}}$ and columns show the same $N_{\mathrm{2}}/N(0)$. The symbol “*” indicates that the model has initial $N(0)=64k$ and “**” indicates that both $N(0)=32k$ and $N(0)=64k$ exist. All other models have initial $N(0)=32k$.  $\frac{m_{\mathrm{2}}}{m_{\mathrm{1}}}\backslash\frac{N_{\mathrm{2}}}{N}(0)$ | 0.005 | 0.01 | 0.03 | 0.05
---|---|---|---|---
2.5 |  | 0.025 | 0.072 | 0.116
5 |  | 0.048 | 0.134 | 0.208
10 | $0.048^{*}$ | $0.092^{**}$ | 0.236 | 0.345
20 |  | $0.168^{**}$ | 0.382 | 0.513
40 |  | $0.288^{**}$ | 0.553 | 0.678

##  4 Isolated star clusters

###  4.1 Energy flux

Figure 1: The evolution of the energy-flux ratios between BH subsystems and
clusters. Upper: $k$ defined in Eq. 5. Middle: the corrected energy-flux
ratio, $k_{\mathrm{1}}=k\psi$ (Eq. 8). Lower: the energy-flux ratio,
$k_{\mathrm{2}}$, by considering both two components inside $R_{\mathrm{h,2}}$
and the relaxation correction factors, $\psi$ and $\psi_{\mathrm{2}}$ (Eq. 9).
The time is in the Hénon unit (NB; Hénon, 1971). All models shown here are in
an isolated environment. Colors represent $m_{\mathrm{2}}/m_{\mathrm{1}}$
(m[value]) and lines indicate $N_{\mathrm{2}}/N(0)$ (N[value]). The models
with initial $64$K particles are shown as $+$ and $\times$ symbols. The two
purple arrows indicate the two stages of evolution: mass segregation and
energy balance.

We first investigate whether the energy-balance theory (Eq. 1) of Breen &
Heggie (2013) applies to star clusters with a large range of
$M_{\mathrm{2}}/M(0)$ in the isolated environment. Since the four quantities
in Eq. 1 can be directly measured in the $N$-body simulations, it is easy to
calculate the energy-flux ratio, $k$. To compute the relaxation time,
$T_{\mathrm{rh}}$, we use the definition from Spitzer (1987),

$T_{\mathrm{r}h}=0.138\frac{N^{1/2}R_{\mathrm{h}}^{3/2}}{\langle
m\rangle^{1/2}G^{1/2}\ln\Lambda},$ (4)

where $N$ is the number of objects, $R_{\mathrm{h}}$ is the half-mass radius,
$\langle m\rangle$ is the average mass, $\ln\Lambda$ is the Coulomb logarithm.
We apply $\ln\Lambda=0.02N$ based on the measurement of Giersz & Heggie
(1996). Instead of the total energy $|E|$ and $|E_{\mathrm{2}}|$, we use the
total kinetic energy $E_{\mathrm{krh}}$ and $E_{\mathrm{krh,2}}$ inside
$R_{\mathrm{h}}$ and $R_{\mathrm{h,2}}$ separately. Thus, Eq. 1 is rewritten
as:

$\frac{E_{\mathrm{krh}}}{T_{\mathrm{rh}}}\approx
k\frac{E_{\mathrm{krh,2}}}{T_{\mathrm{rh,2}}}.$ (5)

This replacement of the energy definition is acceptable when both components
are in virial equilibrium. In Section 4.3, we show that this assumption is
valid for most models. The result of the energy-flux ratio is shown in the
upper panel of Fig. 1. If the balanced evolution has been achieved, $k$ would
keep constant after the mass segregation of BHs. Indeed, the result shows a
two-stage evolution. At the beginning, $k$ decreases due to the effect of mass
segregation. When the balance is established, most models have a constant $k$
until the end of simulations as expected. However, a few models with low
$M_{\mathrm{2}}/M(0)$ (e.g. m5-N0.01 and m10-N0.005) show a growth of $k$. We
name them as “k-growing models”. Besides, we also see that the absolute values
of $k$ at the balanced stage have a large variation depending on
$N_{\mathrm{2}}/N(0)$ and $m_{\mathrm{2}}/m_{\mathrm{1}}$.

To understand the diverged behavior of $k$, we need to consider the two
assumptions involved in Eq. 1, 4 and 5. Firstly, the relaxation time formula
(Eq. 4) is initially derived for a one-component cluster. We assume that by
using the average mass $\langle m\rangle$, $T_{\mathrm{rh}}$ can represent the
average relaxation time for multi-component cases. However, Spitzer & Hart
(1971) showed that this average is not appropriate. It is necessary to
introduce a correction factor, $\psi$, to properly calculate multi-component
relaxation time (priv. comm. M. Gieles and Antonini & Gieles, 2019),

$T_{\mathrm{rh,p}}=\frac{T_{\mathrm{rh}}}{\psi}.$ (6)

This is because the average diffusion coefficient $\langle(\Delta
v_{||})^{2}\rangle$ is not equal to the sum of the coefficients of individual
components (Eq. 16 and 17 in Spitzer & Hart, 1971). This results in the
correction factor for the multi-component relaxation time, $\psi$, defined by:

$\psi=\frac{\sum_{k}n_{\mathrm{k}}m_{\mathrm{k}}^{2}/v_{\mathrm{k}}}{\langle
n\rangle\langle m\rangle^{2}/\langle v\rangle},$ (7)

where $n_{\mathrm{k}}$, $m_{\mathrm{k}}$ and $v_{\mathrm{k}}$ are the number
density, the mass of one object and the mean velocity of the component $k$.
$\langle n\rangle$, $\langle m\rangle$ and $\langle v\rangle$ represent the
average values of all components respectively. Here the definition of $\psi$
is different from the original version (Eq. 24) in Spitzer & Hart (1971),
where the equipartition of energy is assumed. Our definition represents a more
general condition. When all quantities are measured within the same
$R_{\mathrm{h}}$, $n_{\mathrm{k}}$ can be replaced by the total number of
individual components inside $R_{\mathrm{h}}$.

$\psi$ can be directly measured in the $N$-body models. By including $\psi$
and using $T_{\mathrm{rh,p}}$ to replace $T_{\mathrm{rh}}$ in Eq. 5 (shown in
the middle panel of Fig. 1),

$\frac{\psi E_{\mathrm{krh}}}{T_{\mathrm{rh}}}\approx
k_{\mathrm{1}}\frac{E_{\mathrm{krh,2}}}{T_{\mathrm{rh,2}}},$ (8)

the new energy-flux ratio, $k_{\mathrm{1}}=k\psi$, has a common value around
$0.3$ for most models. However, k-growing models still exist.

To explain the growth of $k_{\mathrm{1}}$, we investigate another underlining
assumption. Breen & Heggie (2013) assumes that the cluster center is mostly
occupied by BHs, thus the right-hand side of Eq. 5 only includes the BH
component. To obtain a more accurate energy flux inside $R_{\mathrm{h,2}}$, we
also include the light stars in the calculation of the kinetic energy,
$E_{\mathrm{krh,2a}}$, and the two-component relaxation time with correction
factor, $T_{\mathrm{rh,2a}}/\psi_{\mathrm{2}}$. Thus, the new relation is

$\frac{\psi E_{\mathrm{krh}}}{T_{\mathrm{rh}}}\approx
k_{\mathrm{2}}\frac{\psi_{\mathrm{2}}E_{\mathrm{krh,2a}}}{T_{\mathrm{rh,2a}}},$
(9)

The result is shown in the lower panel of Fig. 1. Now the new energy-flux
ratio has a common stable value, $k_{\mathrm{2}}=0.1-0.2$, for most cases
except the models with $M_{\mathrm{2}}/M(0)>0.35$. When
$M_{\mathrm{2}}/M(0)=1$, where the system become one-component, based on the
definition of Eq. 9, the left and right sides are equivalent, i.e.,
$E_{\mathrm{krh}}\equiv E_{\mathrm{krh,2a}}$, $T_{\mathrm{rh}}/\psi\equiv
T_{\mathrm{rh,2a}}/\psi_{\mathrm{2}}$ and $k_{\mathrm{2}}=1.0$. Thus, it is
expected that there exists a transition region of
$M_{\mathrm{2}}/M(0)\in(\mathcal{R}_{\mathrm{tr}},1.0)$ where $k_{\mathrm{2}}$
is between $0.2$ and $1.0$. The result in Fig. 1 suggests that
$\mathcal{R}_{\mathrm{tr}}\approx 0.35$. In the next section, we show that
this transition criterion is linked to the ratio between the mean velocities
of the two components.

Figure 2: The number of heavy (BH) components, $N_{\mathrm{2}}$, as a function
of time. The plotting style is the same as in Fig. 1.

By including both components in the energy flux calculation, the growth trend
of $k$ and $k_{\mathrm{1}}$ is absent in $k_{\mathrm{2}}$. This suggests that
the energy flux from light components in the cluster center becomes more and
more important in the k-growing models, thus $k$ ($k_{\mathrm{1}}$) increases
as BHs become less dominant. When we consider the escape of BHs, such
behaviour can be naturally explained. In Fig. 2, we show the evolution of the
BH number, $N_{\mathrm{2}}$, for all models. The k-growing models have a small
$N_{\mathrm{2}}$ initially ($320$). Once a large fraction of BHs have escaped,
the interaction between BHs and light stars in the center becomes important.
On the other hand, Fig. 1 shows that the model m10-N0.01 with initially 32k
objects (light blue solid curves) has a growth of $k$ while N64-m10-N0.01
(light blue “$\times$” symbol) has a stable $k$. This also suggests that the
decreasing number of BHs is the reason for the growth trend of $k$
($k_{\mathrm{1}}$).

###  4.2 Correction factor, $\psi$

It is worth checking the values of $\psi$ and $\psi_{\mathrm{2}}$, which
indicate how large the difference of $T_{\mathrm{rh,p}}$ and $T_{\mathrm{rh}}$
can be. Fig. 3 shows the average values of $\psi$ and $\psi_{\mathrm{2}}$
depending on the average mass fraction of BHs, $\langle
M_{\mathrm{2}}/M\rangle$, during the balanced evolution. Notice all quantities
for calculating $\psi$ and $\psi_{\mathrm{2}}$ (Eq. 7) are measured within
$R_{\mathrm{h}}$ and $R_{\mathrm{h,2}}$ separately.

The variation of $\psi$ and $\psi_{\mathrm{2}}$ is large. For the m2.5 models,
where the heavy components represent NS like objects, $\psi$ and
$\psi_{\mathrm{2}}$ are close to unity. When $m_{\mathrm{2}}/m_{\mathrm{1}}$
increases, $\psi$ and $\psi_{\mathrm{2}}$ also increases significantly. There
is no clear trend how $\psi$ and $\psi_{\mathrm{2}}$ depends on $\langle
M_{\mathrm{2}}/M\rangle$ and $m_{\mathrm{2}}/m_{\mathrm{1}}$. This complexity
comes from the additional dependence on the mean velocities of both components
(Eq. 7). When $m_{\mathrm{2}}/m_{\mathrm{1}}\leq 10$, $\psi$ is positively
correlated with $\langle M_{\mathrm{2}}/M\rangle$, while when
$m_{\mathrm{2}}/m_{\mathrm{1}}\geq 20$, there is an anti-correlation. For
$\psi_{\mathrm{2}}$, except the m2.5 models, all other cases show anti-
correlation.

Figure 3: The average correction factor, $\psi$ for $T_{\mathrm{rh}}$ (upper panel) and $\psi_{\mathrm{2}}$ for $T_{\mathrm{rh,2}}$ (lower panel), measured during the balanced evolution. The $x$-axis is the average total mass fraction of BH components, $\langle M_{\mathrm{2}}/M\rangle$, within $R_{\mathrm{h}}$. Colors represent $m_{\mathrm{2}}/m_{\mathrm{1}}$.  Table 3: The Spitzer mass stratification instability criterion $\chi$ for all models at the initial state. $\chi<0.16$ indicates the Spitzer stable case. Thus except the m2.5-N0.1 model, all other models are Spitzer unstable. $\frac{m_{\mathrm{2}}}{m_{\mathrm{1}}}\backslash\frac{N_{\mathrm{2}}}{N}(0)$ | 0.005 | 0.01 | 0.03 | 0.05
---|---|---|---|---
2.5 |  | 0.100 | 0.306 | 0.520
5 |  | 0.565 | 1.729 | 2.942
10 | 1.589 | 3.194 | 9.780 | 16.644
20 |  | 18.069 | 55.325 | 94.150
40 |  | 102.215 | 312.968 | 532.594

The relation between the mean velocities of the two components indicate the
degree of energy equipartition. Whether the two-component system can reach
complete energy equipartition,
$m_{\mathrm{1}}v_{\mathrm{1}}^{2}=m_{\mathrm{2}}v_{\mathrm{2}}^{2}$, depends
on $M_{\mathrm{2}}/M_{\mathrm{1}}$ and $m_{\mathrm{2}}/m_{\mathrm{1}}$. The
criterion is provided by Spitzer (1987):

$\chi=M_{\mathrm{2}}/M_{\mathrm{1}}(m_{\mathrm{2}}/m_{\mathrm{1}})^{3/2}<\beta,$
(10)

where $\beta=0.16$ based on Spitzer’s estimation, but may be different
depending on the structure of the system. Table 3 show $\chi$ for all models
at the initial state. Except the model m2.5-N0.01, all other models have
$\chi>0.16$, i.e.,
$m_{\mathrm{1}}v_{\mathrm{1}}^{2}<m_{\mathrm{2}}v_{\mathrm{2}}^{2}$ during the
balanced evolution. Thus, without knowing the mean velocities,
$v_{\mathrm{1}}$ and $v_{\mathrm{2}}$, it is difficult to estimate $\psi$ only
based on $M_{\mathrm{2}}/M$ and $m_{\mathrm{2}}/m_{\mathrm{1}}$. However, we
can first estimate $\psi$ in two extreme cases. If
$v_{\mathrm{2}}=v_{\mathrm{1}}$, which indicates both components have the same
mean velocity,

$\psi=\frac{n_{\mathrm{1}}m_{\mathrm{1}}^{2}+n_{\mathrm{2}}m_{\mathrm{2}}^{2}}{\langle
n\rangle\langle m\rangle^{2}}.$ (11)

In the case of full energy equipartition,
$m_{\mathrm{1}}v_{\mathrm{1}}^{2}=m_{\mathrm{2}}v_{\mathrm{2}}^{2}$,

$\psi=\frac{n_{\mathrm{1}}m_{\mathrm{1}}^{2.5}+n_{\mathrm{2}}m_{\mathrm{2}}^{2.5}}{\langle
n\rangle\langle m\rangle^{2.5}}.$ (12)

Thus, we can use a power index $\gamma$ to represent the influences from the
ratio of mean velocities. $\psi$ can be written as

$\psi=\frac{n_{\mathrm{1}}m_{\mathrm{1}}^{\gamma}+n_{\mathrm{2}}m_{\mathrm{2}}^{\gamma}}{\langle
n\rangle\langle m\rangle^{\gamma}}=\frac{\langle m^{\gamma}\rangle}{\langle
m\rangle^{\gamma}}.$ (13)

Figure 4: The power index, $\gamma$ for $\psi$ (upper panel) and
$\gamma_{\mathrm{2}}$ for $\psi_{\mathrm{2}}$ (lower panel), from Eq. 13. The
plotting style is the same as in Fig. 3.

From the $N$-body models, we measure the $\gamma$ factor for both $\psi$ and
$\psi_{\mathrm{2}}$. The result is shown in Fig. 4. The $\gamma$ of $\psi$
depends on $\langle M_{\mathrm{2}}/M\rangle$ but is not very sensitive to
$m_{\mathrm{2}}/m_{\mathrm{1}}$. When $\langle M_{\mathrm{2}}/M\rangle$
increases from $0.0$ to $\mathcal{R}_{\mathrm{tr}}$, $\gamma$ decreases from
$2.26$ to $2.0$ and the system departs from energy equipartition. Then
$\gamma$ is stable around $2.0$ and never reach below $2.0$ for $\langle
M_{\mathrm{2}}/M\rangle>\mathcal{R}_{\mathrm{tr}}$. This phenomenon can be
understood when the escape of light stars is considered. The mean escape
velocity of the system can be estimated as

$\langle v_{\mathrm{esc}}\rangle\approx 2\sqrt{|\phi|},$ (14)

where $\phi$ is the total potential of the system, and the mean velocity of
the system in virial equilibrium is

$\langle v\rangle\approx\sqrt{|\phi|}.$ (15)

Thus

$\langle v_{\mathrm{esc}}\rangle\approx 2\langle v\rangle.$ (16)

When $\langle M_{\mathrm{2}}/M\rangle>\mathcal{R}_{\mathrm{tr}}$, heavy
components dominate the potential, energy equipartition leads to a high mean
velocity of light stars. Averagely, light stars with velocities larger than
$\langle v_{\mathrm{esc}}\rangle$ escape the cluster, thus
$v_{\mathrm{1}}<\langle v_{\mathrm{esc}}\rangle$ should be always satisfied.
On the other hand, $v_{\mathrm{2}}$ becomes closer to $\langle v\rangle$ when
$\langle M_{\mathrm{2}}/M\rangle$ increases. Eventually,
$v_{\mathrm{1}}\approx 2v_{\mathrm{2}}$ in the extreme case, which suggests
that $v_{\mathrm{2}}/v_{\mathrm{1}}$ is independent of
$m_{\mathrm{2}}/m_{\mathrm{1}}$ and $\gamma\approx 2.0$. The transition of
$\gamma$ from decrease to flat appears around $\mathcal{R}_{\mathrm{tr}}$,
suggesting that the departure of $k_{\mathrm{2}}$ from $0.1-0.2$ shown in Fig.
1 is linked to the enhanced escape of light stars.

The behavior of $\psi_{\mathrm{2}}$ is similar to $\psi$, but the degree of
energy equipartition inside $R_{\mathrm{h,2}}$ is higher, i.e. the $\gamma$ of
$\psi_{\mathrm{2}}$ ($\gamma_{\mathrm{2}}$) is larger as shown in Fig. 4.
$\gamma_{\mathrm{2}}$ also becomes close to $2.0$ when $\langle
M_{\mathrm{2}}/M\rangle>0.6$. In the center region, the relaxation and mass
segregation times are shorter and the escape velocity (limit for the maximum
mean velocity of light stars) is higher compared to the case in the halo. Thus
it is reasonable to see a higher degree of energy equipartition.
$\gamma_{\mathrm{2}}$ shows a stronger dependence on
$m_{\mathrm{2}}/m_{\mathrm{1}}$, this explains a larger variation of
$\psi_{\mathrm{2}}$ depending on $m_{\mathrm{2}}/m_{\mathrm{1}}$.

Figure 5: The ratio between the light-component relaxation time,
$T_{\mathrm{rh,1}}$, and the corrected two-component relaxation time,
$T_{\mathrm{rh,p}}$. The plotting style is the same as in Fig. 3.

Observations of star clusters cannot directly measure the two-component
relaxation time, $T_{\mathrm{rh,p}}$, since BHs are invisible. Thus we
calculate the ratio between light-component relaxation time,
$T_{\mathrm{rh,1}}$, and $T_{\mathrm{rh,p}}$ in Fig. 5. This ratio
significantly increases with larger $M_{\mathrm{2}}/M$ and varies in a range
of $1-70$. For $M_{\mathrm{2}}/M>0.4$, BHs significantly dominate the system,
thus the ratio can be above $20$. In other cases, the ratio is below $10$.
This result indicates that by only measuring the relaxation time of stars, we
would significantly overestimate the two-component relaxation time if
$M_{\mathrm{2}}/M$ is large.

###  4.3 Half-mass radius

Figure 6: The evolution of $k^{\prime}$ defined by Eq. 17 (upper panel). The
coefficients with correction $k^{\prime}_{\mathrm{1}}=k^{\prime}\psi$ are
shown in the lower panel. The plotting style is the same as in Fig. 1.

The importance of $\psi$ should also reflect in the relation of
$R_{\mathrm{h,2}}/R_{\mathrm{h}}$ shown in the Eq. 4 of Breen & Heggie (2013):

$\left(\frac{R_{\mathrm{h,2}}}{R_{\mathrm{h}}}\right)^{5/2}\approx
k^{\prime}\left(\frac{M_{\mathrm{2}}}{M}\right)^{3/2}\frac{m_{\mathrm{2}}}{\langle
m\rangle}\frac{\ln\Lambda_{\mathrm{2}}}{\ln\Lambda}.$ (17)

This is derived based on Eq. 1. When the two components are in virial
equilibrium, $k^{\prime}=k$.

In the upper panel of Fig. 6, we compare the coefficient $k^{\prime}$ of
different models defined by Eq. 17. As expected, the differences of $k$ shown
in Fig. 1 also exist for $k^{\prime}$. For the models with flat $k^{\prime}$
after the balanced evolution start, $k\approx k^{\prime}$, suggesting that the
two components are in virial equilibrium. The models with increasing
$k^{\prime}$ have $k^{\prime}<k$ after about $8000$ NB time unit, which
suggests a departure from virial equilibrium.

In the lower panel, we apply the correction $\psi$ like the case of Eq. 8, the
differences are reduced, similar to the middle panel of Fig. 1. Thus, by using
Eq. 13, a new relation for the half-mass-radius ratio can be written as

$\begin{split}\left(\frac{R_{\mathrm{h,2}}}{R_{\mathrm{h}}}\right)^{5/2}\approx&k^{\prime}_{\mathrm{1}}\left(\frac{M_{\mathrm{2}}}{M}\right)^{1/2}\left(\frac{m_{\mathrm{2}}}{\langle
m\rangle}\right)^{2-\gamma}\\\
&\frac{1}{1+\frac{M_{\mathrm{1}}}{M_{\mathrm{2}}}\left(\frac{m_{\mathrm{1}}}{m_{\mathrm{2}}}\right)^{\gamma-1}}\frac{\ln\Lambda_{\mathrm{2}}}{\ln\Lambda},\end{split}$
(18)

To obtain this relation, the number densities $n_{\mathrm{1}}$,
$n_{\mathrm{2}}$ and $\langle n\rangle$ in Eq. 13 are replaced by
$M_{\mathrm{1}}/m_{\mathrm{1}}$, $M_{\mathrm{2}}/m_{\mathrm{2}}$ and
$M/\langle m\rangle$ separately (measured within $R_{\mathrm{h}}$). Since the
corrections of relaxation time and energy flux due to the existence of light
stars inside the region of BH subsystems (Eq. 9) are not included here, this
relation only provides a better but not exactly accurate description of
$R_{\mathrm{h,2}}/R_{\mathrm{h}}$ compared to Eq. 17. It is valid when
$M_{\mathrm{2}}/M<\mathcal{R}_{\mathrm{tr}}$ and $N_{\mathrm{2}}>300$ (a rough
boundary).

###  4.4 Mass loss

The dissolution time of star clusters depends on the relaxation time. In this
section, we study whether the corrected relaxation time, $T_{\mathrm{rh,p}}$,
is indeed consistent with the dissolution time of clusters. In Fig. 7 and 8,
we compare the evolution of the remaining masses of different components
normalized to their initial values, $M_{\mathrm{2}}(t)/M_{\mathrm{2}}(0)$
(BHs), $M_{\mathrm{1}}(t)/M_{\mathrm{1}}(0)$ (stars) and $M(t)/M(0)$ (both),
based on different time units. In each sub-plot, the time is scaled in the way
of an integration:

$T[T_{\mathrm{scale}}]=\int^{t}_{t_{\mathrm{0}}}\frac{dt}{T_{\mathrm{scale}}(t)},$
(19)

where the choice of $T_{\mathrm{scale}}$ is shown in the $x$-axis labels of
sub-plots. The starting time, $t_{\mathrm{0}}=10~{}T[T_{\mathrm{ms}}]$, is the
roughly measured time of BH core collapse (or the finishing time of BH mass
segregation; except for the m2.5 models). $T_{\mathrm{ms}}$ is the mass-
segregation timescale of BHs (Spitzer, 1987):

$T_{\mathrm{ms}}=\frac{m_{\mathrm{1}}}{m_{\mathrm{2}}}T_{\mathrm{rh,1}}.$ (20)

The normalization factors of the remaining masses, $M_{\mathrm{1}}(0)$,
$M_{\mathrm{2}}(0)$ and $M(0)$, are also determined at $t_{\mathrm{0}}$, thus
all models have the y-axis value of $1.0$ when $t=t_{\mathrm{0}}$. This choice
of $t_{\mathrm{0}}$ makes a better comparison of the mass loss among different
models, because after core collapse, mass loss becomes significant, especially
for BHs. If the choice of $T_{\mathrm{scale}}$ provides a proper description
of the mass-loss timescale, the curves (normalized remaining masses) from
different models are expected to overlap with each other.

####  4.4.1 BHs

Figure 7: The evolution of the remaining mass of BHs normalized to their
initial values. From left to right, the four columns show the time in the
units of $T_{\mathrm{rh,2}}$ (the half-mass relaxation time of BH subsystems),
$T_{\mathrm{rh}}$ (the average half-mass relaxation time of clusters),
$T_{\mathrm{ms}}=(m_{\mathrm{1}}/m_{\mathrm{2}})T_{\mathrm{rh,1}}$ (the BH
mass-segregation timescale) and $T_{\mathrm{rh,p}}=T_{\mathrm{rh}}/\psi$ (the
corrected two-component half-mass relaxation time) separately (see Eq. 19).
For better comparison, $t=0$ is defined at $10~{}T[T_{\mathrm{ms}}]$ after the
starting time of simulations. The dashed line indicates
$M_{\mathrm{2}}/M_{\mathrm{2}}(0)=0.5$. The plotting style is the same as in
Fig. 1.

Breen & Heggie (2013) found that the ejections via few-body interaction should
be the major reason for the escape of BHs. The relaxation driven evaporation
mechanism is negligible because most BHs are trapped in the deep potential of
the cluster. They assume that encounters with one hard BH binary in the
cluster core release a certain amount of energy before the binary escape the
cluster, and this energy generation rate is controlled by the global energy
requirement. On the other hand, each BH binary ejected a fixed number ($5-6$;
including the binary itself) of BHs (Goodman, 1984; Heggie & Hut, 2003).
Therefore, the escape rate of BHs depends on the global half-mass relaxation
time, $T_{\mathrm{rh}}$, and the total mass, $M$, and is independent on
$T_{\mathrm{rh,2}}$, as shown in Eq. 9 of Breen & Heggie (2013)

$\dot{M}_{\mathrm{2}}\sim-M/T_{\mathrm{rh}}.$ (21)

This relation is derived for the case of $M_{\mathrm{2}}/M\approx 0.01$, where
the central potential is dominated by the light components. Consider $\psi>1$
when $M_{\mathrm{2}}/M$ is large, it is expected that $T_{\mathrm{rh,p}}$
should replace $T_{\mathrm{rh}}$, thus

$\dot{M}_{\mathrm{2}}\sim-MT_{\mathrm{rh,p}}=-M\psi/T_{\mathrm{rh}}.$ (22)

Fig. 7 show how $M_{\mathrm{2}}(t)/M_{\mathrm{2}}(0)$ depends on different
scaled time. Four types of $T_{\mathrm{scale}}$ are compared:
$T_{\mathrm{rh,2}}$, $T_{\mathrm{rh}}$, $T_{\mathrm{ms}}$ and
$T_{\mathrm{rh,p}}$. The loss of half BHs takes $100-1000$ $T_{\mathrm{rh,2}}$
or a few $T_{\mathrm{rh}}$, suggesting that the survival timescale of BHs is
much longer than $T_{\mathrm{rh,2}}$ and is close to $T_{\mathrm{rh}}$ or
$T_{\mathrm{rh,p}}$. This is consistent with the expectation of Eq. 21 or 22.
However, when time is in the units of $T_{\mathrm{rh}}$ or
$T_{\mathrm{rh,p}}$, the curves $M_{\mathrm{2}}(t)/M_{\mathrm{2}}(0)$ show a
significant variation of slopes depending on $N_{\mathrm{2}}/N$ and
$m_{\mathrm{2}}/\langle m\rangle$, This indicates that there is a missing
coefficient that depends on the properties of BH subsystems in Eq. 21 and 22.
Instead, the 3rd panel in Fig. 7 shows a well consistent evolution (same
slopes) of $M_{\mathrm{2}}(t)/M_{\mathrm{2}}(0)$ when time is in the unit of
$T_{\mathrm{ms}}$. The mass loss of BH subsystems can be approximately
described by:

$\begin{split}M_{\mathrm{2}}(t)/M_{\mathrm{2}}(t_{\mathrm{0}})&\approx
0.0125\int^{t}_{t_{\mathrm{0}}}\frac{dt}{T_{\mathrm{ms}}}\\\ \end{split}$ (23)

Now we try to figure out why the expectation of Eq. 22 is not exactly
reflected in the numerical models. The idea that one binary eject a fixed
number of stars is obtained from the analysis of one-component systems. In the
two-component case, a BH binary ejects both BHs and stars. In Section 4.1, we
have shown that the light stars are necessary to be included in the
calculation of the energy flux in the cluster center. Thus, the ejection of
stars cannot be completely ignored. Then, the number of ejected BHs per BH
binary depends on the number densities of two components in the core, which
rely on $N_{\mathrm{2}}/N$ and $m_{\mathrm{2}}/\langle m\rangle$. This
dependence is hidden in the missing coefficient of Eq. 22, which causes the
diverged behavior of models shown in Fig. 7.

To understand why mass segregation plays the dominant role, we start from a
different point of view: the idea of “loss cone” from the the studies of
supermassive black hole binaries (SMBHBs) in the galactic nuclei (e.g.
Lightman & Shapiro, 1977). The loss cone is the region of phase space where
stars are on orbits intersecting the orbit of SMBHBs. These stars are quickly
ejected and the loss cone is depleted. A similar situation can happen in the
case of a hard BH binary in the cluster center. To support the energy
requirement, we assume that the core contains at least one hard BH binary in
most of the time after core collapse. If initially the loss cone of the hard
BH binary is empty, it can be refilled by BHs suffering orbital changes via
perturbations. These BHs are quickly ejected after an strong interaction with
the binary. In two-component systems, mass segregation is more efficient than
two-body relaxation to refill the loss cone. If a fixed fraction of the
ejected BHs exceed the escape velocity of the cluster, it is expected to see
that the mass-loss rate of BHs depends on $T_{\mathrm{ms}}$.

####  4.4.2 Light stars

Figure 8: The evolution of the remaining masses of light stars (upper) and
both (lower) normalized to their initial values separately. The two columns
show the time in the units of $T_{\mathrm{rh}}$ and $T_{\mathrm{rh,p}}$
separately. The plotting style is the same as in Fig. 7.

In the isolated models, the escape of light stars is driven by the two-body
relaxation evaporation and ejections after strong interactions with binaries.
The former is considered as the major channel, thus the mass-loss rate of
light stars depends on the relaxation time. We expect that $T_{\mathrm{rh,p}}$
should provide a better match to the mass-loss rate of light stars. We compare
the evolution of the normalized light-component mass,
$M_{\mathrm{1}}(t)/M_{\mathrm{1}}(0)$, and the normalized total mass,
$M(t)/M(0)$, with different scaled times in Fig. 8 (similar like the analysis
in Section 4.4.1). The two columns show the times in $T_{\mathrm{rh}}$ and
$T_{\mathrm{rh,p}}$ separately. Although a large divergence appears in both
cases, $T[T_{\mathrm{rh,p}}]$ indeed provides a better result.

There is also a clear difference for models with $M_{\mathrm{2}}/M$ below and
above $\mathcal{R}_{\mathrm{tr}}$. The models of
$M_{\mathrm{2}}/M<\mathcal{R}_{\mathrm{tr}}$ shows a much slower dissolution.
In Section 4.2, we show that when
$M_{\mathrm{2}}/M>\mathcal{R}_{\mathrm{tr}}$, the energy equipartition drives
the mean velocity of stars increasing close to the mean escape velocity of the
system. Therefore the escape of light stars is expected to be accelerated when
BHs dominate the system. The speed up of cluster dissolution due to the BH
subsystem is also observed in the Monte-Carlo simulations of GCs (Giersz, et
al., 2019).

Figure 9: The ratio between the minimum ejection velocity of a star
($v_{\mathrm{ej,m}}$) from an encounter with a BH binary and the central
escape velocity of clusters ($v_{\mathrm{esc,c}}$). The time is in the unit of
$T_{\mathrm{rh,p}}$. The plotting style is the same as in Fig. 7.

On the other hand, we can estimate the ejection velocity of a light star,
$v_{\mathrm{ej}}$, after an strong interaction with a BH binary in the cluster
core by (Miller & Lauburg, 2009; Antonini, Gieles, & Gualandris, 2019):

$v_{\mathrm{ej}}\approx\left[0.2G\frac{m_{\mathrm{1}}m_{\mathrm{2}}}{2(m_{\mathrm{1}}+2m_{\mathrm{2}})a}\right]^{1/2}.$
(24)

where $a$ is the semi-major axis of the binary. By using $a$ at the half-soft
boundary (Heggie-Hills law; Heggie, 1975; Hills, 1975) in BH subsystems,

$a_{\mathrm{hs}}=\frac{Gm_{\mathrm{2}}^{2}}{\langle m\rangle
v_{\mathrm{2}}^{2}},$ (25)

the lower boundary of $v_{\mathrm{ej,m}}$ can be estimated. Considering the
mean escape velocity for the Plummer model at the cluster center,

$v_{\mathrm{esc,c}}=\left(2\times
1.305G\frac{M}{R_{\mathrm{h}}}\right)^{1/2},$ (26)

we now calculate the ratio, $v_{\mathrm{ej,m}}/v_{\mathrm{esc,c}}$, as shown
in Fig. 9. Interestingly, the models with
$M_{\mathrm{2}}/M>\mathcal{R}_{\mathrm{tr}}$ have
$v_{\mathrm{ej,m}}/v_{\mathrm{esc,c}}>1.0$. Thus, it is easier for a star to
escape the cluster by one strong binary-single interaction. This may be
another important reason for the fast dissolution.

##  5 Tidally filling star clusters

###  5.1 Energy flux

Figure 10: The evolution of the energy-flux ratio for tidally filling models,
similar to Fig. 1. The time is in the unit of the dissolution time,
$T_{\mathrm{dis,p}}(x=1)$, defined in Eq. 29.

When star clusters are in the galactic tidal field, the tidal stripping
becomes an important mechanism to remove stars and BHs from the cluster. We
first investigate how the tidal field influences the balance of energy flux.
Fig. 10 shows the evolution of energy flux ratio (Eq. 9), similar to Fig. 1.
Instead of keeping a constant value as in the isolated star clusters,
$k_{\mathrm{2}}$ increases for all models after mass segregation. For most
cases, $k_{\mathrm{2}}$ first decreases to $0.1\sim 0.2$ and then increases to
about $1.0$ at the end. The increasing $k_{\mathrm{2}}$ is not necessary to
suggest that the energy balance is lost. Instead, it indicates that the two-
component system is smoothly transferred to an one-component system (we also
discuss this point in Section 5.4). When a significant fraction of light stars
are removed via tidal stripping, $R_{\mathrm{h,2}}$ increases towards
$R_{\mathrm{h}}$ and the energy flux measured within $R_{\mathrm{h,2}}$ does
not represent the energy generation at the cluster center any more.

###  5.2 Mass loss of BHs

Figure 11: The normalized remaining mass of BHs for tidally filling models
depending on the time in the unit of $T_{\mathrm{ms}}$. The plotting style is
the same as in Fig. 7.

Before the halo of light stars are evaporated, BHs are inside the deep
potential. During this period, the escape rate of BHs is expected to be
independent of tidal fields and be similar to the case of isolated clusters,
i.e., can be described by Eq. 23. In Fig. 11, we find that
$M_{\mathrm{2}}(t)/M_{\mathrm{2}}(0)$ for tidally filling star clusters is
still consistent with $T[T_{\mathrm{ms}}]$, but the models with high
$M_{\mathrm{2}}/M$ lose BHs faster. This is because the high
$M_{\mathrm{2}}/M$ causes faster expansion of $R_{\mathrm{h,1}}$, which
indicates that light stars are easier to escape. Eventually when a large
fraction of light stars escape, the tidal stripping of BHs becomes more and
more important.

###  5.3 Dissolution of clusters

Figure 12: The evolution of total mass normalized to the initial values for
tidally filling models. The five panels show the time in the units of
dissolution times defined as: $T_{\mathrm{dis}}(x=1)$,
$T_{\mathrm{dis,p}}(x=1)$, $T_{\mathrm{dis,p}}(x=0.75)$,
$T_{\mathrm{dis,p1}}(x=0.75)$ and $T_{\mathrm{dis,p}}(x=1;init)$ separately.
$T_{\mathrm{dis}}$ (Eq. 2) use the average relaxation time, $T_{\mathrm{rh}}$,
while $T_{\mathrm{dis,p}}$ and $T_{\mathrm{dis,p1}}$ use the corrected value,
$T_{\mathrm{rh,p}}$. $x=1$ indicates no back-scatter effect of potential
escapers, while $x=0.75$ indicates the opposite case. The last panel shows the
time in the unit of the initial value of $T_{\mathrm{dis,p}}(x=1)$. The
plotting style is the same as in Fig. 7.

The dissolution time of tidal filling star clusters with one-component or
Kroupa/Chabrier IMFs are well described by Eq. 2 (Baumgardt, 2001; Baumgardt &
Makino, 2003). However, the maximum mass of zero-age main sequence stars in
their $N$-body models is $15M_{\odot}$. Thus BH populations are excluded in
their analysis. For the Kroupa/Chabrier IMFs, the final population of BHs is a
small fraction of the total mass (see Table 1), the difference between
$T_{\mathrm{rh,p}}$ and $T_{\mathrm{rh}}$ is small, i.e., the global evolution
of the cluster is dominated by the light component. In our models that
represent clusters with top-heavy IMFs, this difference cannot be ignored. By
using $T_{\mathrm{rh,p}}$ instead of $T_{\mathrm{rh}}$ in Eq. 2, the corrected
dissolution time is

$T_{\mathrm{dis,p}}\sim
T_{\mathrm{rh,p}}^{x}T_{\mathrm{cr}}^{1-x}(R_{\mathrm{t}}/R_{\mathrm{h}})^{3/2}.$
(27)

But there is one complexity in the two-component models when calculating the
radius ratio, $R_{\mathrm{t}}/R_{\mathrm{h}}$, and the crossing time,
$T_{\mathrm{cr}}$. The latter represents the timescale for a star escaping the
system after gaining enough energy to become a potential escaper. However, the
heavy and light components have different half-mass radii and crossing times
in the balanced-evolution phase. Especially when $M_{\mathrm{2}}/M$ is large,
the half-mass radius of stars, $R_{\mathrm{h,1}}$, can be much larger than
$R_{\mathrm{h}}$ of the system. Thus, $T_{\mathrm{cr}}$ calculated by using
the averaged property of the system may not correctly represent the escape
time of light stars. To investigate this, we define another dissolution time
by using $R_{\mathrm{h,1}}$ and $T_{\mathrm{cr,1}}$ (the crossing time
measured at $R_{\mathrm{h,1}}$):

$T_{\mathrm{dis,p1}}\sim
T_{\mathrm{rh,p}}^{x}T_{\mathrm{cr,1}}^{1-x}(R_{\mathrm{t}}/R_{\mathrm{h,1}})^{3/2}.$
(28)

In Fig. 12, we compare the evolution of the total mass, $M(t)/M(0)$, depending
on the time in the units of $T_{\mathrm{dis}}$, $T_{\mathrm{dis,p}}$ and
$T_{\mathrm{dis,p1}}$ with $x=1$ and $x=0.75$. The dissolution time with $x=1$
excludes the back-scatter effect of potential escapers while $x=0.75$
represents the opposite case. The result shown in Fig. 12 indicates that
$T_{\mathrm{dis,p}}(x=1)$ is the best timescale for describing the cluster
dissolution (see the 2nd panel). By including the back-scatter effect, a large
divergence appears among different models, as shown in the 3rd panel. But once
the crossing time and half-mass radius use the values measured at
$R_{\mathrm{h,1}}$ (the 4th panel), the behaviour of different models becomes
much more consistent. This indicates the importance to have a proper
definition of the crossing time and the radius scaling factor in the formula
of the dissolution time. If the time is in the unit of the initial value of
$T_{\mathrm{dis,p}}(x=1)$, shown in the 5th panel, the decrease of $M(t)/M(0)$
becomes more linear, but the divergence of different models also become
larger. This is because the dissolution time evolves significantly due to an
expansion of the system and a large change of $M_{\mathrm{2}}/M$. By using the
result of $T_{\mathrm{dis,p}}(x=1)$, we obtain

$M(t)/M(t_{\mathrm{0}})\approx
0.5\int^{t}_{t_{\mathrm{0}}}\frac{dt}{T_{\mathrm{dis,p}}(x=1)}.$ (29)

This is valid when the cluster is initially tidally filling.

###  5.4 The evolution of $M_{\mathrm{2}}/M$

Eq. 23 and 29 show that the mass loss of BH subsystems and clusters have
different timescales. The BHs escape on the timescale of mass segregation
(Fig. 7 and 11), but the mass loss of light stars depends on the relaxation
time (Fig. 8 and 12). On the other hand, the dissolution of light components
relies on tidal fields, while the escape of BHs does not when they are deeply
trapped in the cluster center. Thus, whether $M_{\mathrm{2}}/M$ increases or
decreases during the long-term evolution is sensitive to the properties of BH
subsystems (IMFs) and tidal fields. Fig. 13 shows that without and with a
tidal field, the evolution of $M_{\mathrm{2}}/M$ is very different.

The isolated models with initial $M_{\mathrm{2}}/M<\mathcal{R}_{\mathrm{tr}}$
have a decreasing trend of $M_{\mathrm{2}}/M$, which indicates that light
stars become the dominant component and only few or no BHs stay in clusters
after a long-term evolution. For other models, $M_{\mathrm{2}}/M$ slowly
increases, eventually the systems evolve to dark clusters without light stars.
As discussed in Section 4.2 and 4.4, when
$M_{\mathrm{2}}/M>\mathcal{R}_{\mathrm{tr}}$, the escape of light stars is
more efficient as a result of energy equipartition. Thus, it is expected that
there is such a transition criterion determining the fate of clusters.

In tidally filling models, due to the tidal stripping of light stars, the
evolution of $M_{\mathrm{2}}/M$ becomes faster. The transition boundary of
$\mathcal{R}_{\mathrm{tr}}$ also becomes lower. $M_{\mathrm{2}}/M$ in the
models with initial $M_{\mathrm{2}}/M(0)>0.07$ increases fast and finally
becomes $1.0$, which means that the systems evolve to BH-dominant dark
clusters. This transition boundary is also discussed in Breen & Heggie (2013),
where they suggested a constant value that $\mathcal{R}_{\mathrm{tr}}\approx
0.11$. The fixed value of $0.11$ is obtained based on the assumption that the
mass loss of light stars and BHs both depend on the relaxation time and the
total mass. Although $0.11$ is not far from the value in our result, the
boundary, $\mathcal{R}_{\mathrm{tr}}$, is probably not a constant but depends
on $m_{\mathrm{2}}/m_{\mathrm{1}}$, $M_{\mathrm{2}}/M$ and
$R_{\mathrm{t}}/R_{\mathrm{h}}$. Fig. 13 shows that the slope of
$M_{\mathrm{2}}/M$ for different models is not the same, which reflects such
dependence. On the other hand, the result in Fig. 13 indicates that when a
star cluster has an IMF with $\alpha_{\mathrm{3}}>-2.0$, where
$M_{\mathrm{2}}/M(0)\geq 0.182$, it can easily evolve to a dark cluster after
about $20~{}T_{\mathrm{rh,p}}$ in a strong tidally filling environment.

Figure 13: The evolution of $M_{\mathrm{2}}/M$ for tidally filling and
isolated models. The plotting style is the same as in Fig. 7.

##  6 Discussion and conclusion

In this work, we have carried out a series of $N$-body simulations of star
clusters with two components (BHs and light stars), which approximates star
clusters with top-heavy IMFs. The simplification of the $N$-body models allows
us to focus on the dynamical effect of BH subsystems in the long-term
evolution of star clusters. We validate the idea of energy-flux balanced
evolution based on the Hénon (1975) principle and Breen & Heggie (2013). The
result (Fig. 1) indicates that to properly calculate the two-body relaxation
time of two-component clusters, a correction factor $\psi$ to the Spitzer
relaxation time for average cluster properties, $T_{\mathrm{rh}}$ (Eq. 4), is
necessary ($T_{\mathrm{rh,p}}=T_{\mathrm{rh}}/\psi$). Our models show that
$\psi=1-10$ for star clusters with the total mass fraction of BHs,
$M_{\mathrm{2}}/M$, in the range of $2.5\%-68\%$. Since $\psi$ depends on the
property of BHs (Eq. 7 and 13), the evolution of star clusters, especially the
dissolution time (Eq. 27 and 28), are sensitive to $M_{\mathrm{2}}/M$ and the
mass ratio, $m_{\mathrm{2}}/m_{\mathrm{1}}$. Thus, star clusters with top-
heavy IMFs tend to dissolve faster. For example, with
$\alpha_{\mathrm{3}}=-1.7$ (see Table 2), the initial $M_{\mathrm{2}}/M\approx
46.6\%$, the corresponding $\psi\approx 10$. It is expected that such clusters
can dissolve $10$ times faster than the estimation based on the relaxation
time for average cluster properties. Moreover, in observations BHs cannot be
directly detected, thus the measurement of relaxation time only counts light
stars. The comparison of $T_{\mathrm{rh,p}}$ and the light-component
relaxation time, $T_{\mathrm{rh,1}}$, shows that $T_{\mathrm{rh,1}}$ can be
$1-70$ times larger than $T_{\mathrm{rh,p}}$. This suggests that if we ignore
BHs in a star cluster, the relaxation time of the system can be significantly
overestimated.

The $\psi$ factor can be simplified as Eq. 13 with a power index $\gamma$. The
value of $\gamma$ represents the degree of mass segregation of two-component
systems. In Fig. 3 and 4, we obtain the average $\psi$ and $\gamma$ depending
on $m_{\mathrm{2}}/m_{\mathrm{1}}$ and $\langle M_{\mathrm{2}}/M\rangle$. The
average is calculated by using data from the start of the balanced evolution
to the end of simulations. For models with small values of $\langle
M_{\mathrm{2}}/M\rangle$ which represent the case of the canonical IMF, the
system has a high degree of mass segregation with $\gamma$ being close to
$2.5$. When $M_{\mathrm{2}}/M>0.35$ (top-heavy IMF), the mean velocities of
two components tend to be close to each other and be independent of
$m_{\mathrm{2}}/m_{\mathrm{1}}$ ($\gamma\approx 2.0$), because the mean
velocity of light stars is limited by the mean escape velocity of the system.

The escape rate of BHs depends on the mass segregation time (Eq. 23, Fig.
8,11) and is not sensitive to the tidal field (when light stars still dominate
the halo of clusters). This is different from the mass-loss rate of light
stars (Fig. 8 and 12). In the isolated cluster, $T_{\mathrm{rh,p}}$ is a
better timescale than $T_{\mathrm{rh}}$ to describe the mass loss of light
stars. However, the dissolution time shows a very diverged behavior between
clusters with low and high $M_{\mathrm{2}}/M$. In the case of high
$M_{\mathrm{2}}/M$, light stars are easier to reach the escape criterion due
to the energy equipartition and strong interactions with BH binaries. In the
tidal filling case, the dissolution time of the system can be well described
by $T_{\mathrm{rh,p}}$ (Eq. 27) without the back-scatter effect of potential
escaper ($x=1$) or with the effect but only if the crossing time and radius
scaling factor is properly chosen (Eq. 28).

Because the mass-loss rate of BHs and light stars depends on different
timescales, the initial $M_{\mathrm{2}}/M$ and the strength of tidal fields
determine the fate of star clusters. Roughly, with a top-heavy IMF that
$\alpha_{\mathrm{3}}<-2.0$ and a weak tidal field, they eventually evolve to
“luminous clusters” where very few BHs exist. With $\alpha_{\mathrm{3}}>-2.0$
or a strong tidal field, “dark clusters” dominated by BHs easily form. For
isolated star clusters, the boundary of $\langle M_{\mathrm{2}}/M\rangle$ for
this transition is around $\mathcal{R}_{\mathrm{tr}}=0.3-0.4$. For tidally
filling clusters, it is about $0.07$, which is the case of the canonical IMF.
The timescale to become a dark cluster is about $20~{}T_{\mathrm{rh,p}}$ in a
strong tidal filling environment. This suggests that if massive star clusters
(e.g. GCs) with top-heavy IMFs have formed at the high redshift and have been
close to the Galactic center, they were easier to be evaporated by tidal
fields or have become “dark clusters” that cannot be observed by optical
telescopes today. On the other hand, the post-core collapse GCs with dense
cores (e.g. $47$ Tuc; Hénault-Brunet, et al., 2019) do not tend to have top-
heavy IMFs (unless most BHs have escaped after supernovae with high-velocity
natal kicks). The future observations of gravitational lensing and
gravitational waves may provide a constraint on the population of dark
clusters.

In this work, we only theoretically focus on the dynamical effect of BH
subsystems for the long-term evolution of two-component, spherical and low-$N$
star clusters. These simple models allow us to well isolate the dynamical
effect of BH subsystems from other physical processes. In a realistic
situation, primordial binaries, stellar evolution, mass spectrum, evolution of
tidal field, density profiles and rotation of clusters complicate the simple
picture discussed here (e.g. Giersz, et al., 2019).

Especially, the mass spectrum of stars and BHs introduces an additional
complexity in determining the proper relaxation time. We assume that the equal
mass of BHs (stars) can reasonably represents the averaged properties of the
objects with a mass spectrum. However, the importance of $\psi$ discussed in
this work suggests that such approximation may not be accurate when the mass
spread is wide. From Fig. 3, we see that when
$m_{\mathrm{2}}/m_{\mathrm{1}}>5$, $\psi$ becomes larger than $1.0$. In an
evolved star cluster, the mass spectrum of light stars can ranges from $<0.08$
(brwon dwarfs) to $\sim 1.0$ M⊙, thus $\psi$ for only light stars may already
be larger than $1.0$. In such case, the calculation of an accurate relaxation
time requires an integration of the diffusion coefficients of different mass
components. However, Baumgardt & Makino (2003) indicated that a star cluster
with a mass spectrum but no BHs can be well described by Eq. 2 without the
correction factor ($\psi$). This suggests that compared to the strong impact
of BH subsystems shown in this work, the effect introduced by the mass
spectrum of light stars is not significant. Thus the two-component system
should provide a reasonable representation of a star cluster with a top-heavy
IMF.

On the other hand, the stellar evolution brings another complexity.
Especially, the stellar evolution models have large uncertainty, how does the
kick velocity distribution of BHs and NSs after supernovae looks like is still
an open question. But this distribution is crucial to determine the remaining
total mass fraction of BHs and NSs. This leads to an uncertainty to map an IMF
to the initial values of $M_{\mathrm{2}}/M$ and
$m_{\mathrm{2}}/m_{\mathrm{1}}$. Besides, $m_{\mathrm{2}}/m_{\mathrm{1}}$ also
depends on the metallicity. Thus, our results only provide an approximate
analysis to show how to connect the shape of IMFs to the dynamical evolution
of star clusters.

The existence of primordial binaries provides an additional energy source and
increases the possibility of binary-binary interactions. It may influence the
core radius of clusters and the energy generation. However, we expect the
global evolution (mass loss and energy balance) may not be significantly
influenced because whatever the source of energy generation looks like, the
energy generation rate should be always guided by the global energy
requirement (Hénon’s principle).

We focus on the long-term evolution of star clusters. In the formation stage
of clusters embedded in the gas cloud, the influence of stellar-wind feedback
from massive OB stars to the gas can be strongly enhanced when IMFs become
top-heavy. Thus, whether the star formation region with a top-heavy IMF can
become a gravitational bound cluster or an association after gas is removed is
an open question. In the follow-up projects, it is necessary to validate our
conclusion by including these multiple physical effects step by step.

## Acknowledgments

L.W. thanks the Alexander von Humboldt Foundation for funding this research.
L.W. also thanks Douglas C. Heggie, Mark Gieles, Mirek Giersz, Anna Lisa Varri
and Sverre Aarseth for very useful discussions and advice for this project.

## References

  * Aarseth (2003) Aarseth S. J., 2003, Gravitational N-Body Simulations, Cambridge University Press
  * Alexander & Gieles (2012) Alexander P. E. R., Gieles M., 2012, MNRAS, 422, 3415
  * Antonini, Gieles, & Gualandris (2019) Antonini F., Gieles M., Gualandris A., 2019, MNRAS, 486, 5008
  * Antonini & Gieles (2019) Antonini F., Gieles M., 2019, arXiv, arXiv:1906.11855
  * Banerjee & Kroupa (2011) Banerjee S., Kroupa P., 2011, ApJL, 741, L12
  * Banerjee, et al. (2019) Banerjee S., Belczynski K., Fryer C. L., Berczik P., Hurley J. R., Spurzem R., Wang L., 2019, arXiv, arXiv:1902.07718
  * Bastian, Covey & Meyer (2010) Bastian N., Covey K. R., Meyer M. R., 2010, ARA&A, 48, 339
  * Baumgardt (2001) Baumgardt H., 2001, MNRAS, 325, 1323
  * Baumgardt & Makino (2003) Baumgardt H., Makino J., 2003, MNRAS, 340, 227
  * Baumgardt & Sollima (2017) Baumgardt H., Sollima S., 2017, MNRAS, 472, 744
  * Belczynski, et al. (2016) Belczynski K., et al., 2016, A&A, 594, A97
  * Breen & Heggie (2012) Breen P. G., Heggie D. C., 2012, MNRAS, 420, 309
  * Breen & Heggie (2013) Breen P. G., Heggie D. C., 2013, MNRAS, 432, 2779
  * Chabrier (2003) Chabrier G., 2003, PASP, 115, 763
  * Chatterjee, Rodriguez, & Rasio (2017) Chatterjee S., Rodriguez C. L., Rasio F. A., 2017, ApJ, 834, 68
  * Fryer, et al. (2012) Fryer C. L., Belczynski K., Wiktorowicz G., Dominik M., Kalogera V., Holz D. E., 2012, ApJ, 749, 91
  * Goodman (1984) Goodman J., 1984, ApJ, 280, 298
  * Gieles, Heggie, & Zhao (2011) Gieles M., Heggie D. C., Zhao H., 2011, MNRAS, 413, 2509
  * Hénault-Brunet, et al. (2019) Hénault-Brunet V., Gieles M., Strader J., Peuten M., Balbinot E., Douglas K. E. K., 2019, MNRAS.tmp, 2599
  * Giersz & Heggie (1996) Giersz M., Heggie D. C., 1996, MNRAS, 279, 1037
  * Giersz, et al. (2019) Giersz M., Askar A., Wang L., Hypki A., Leveque A., Spurzem R., 2019, MNRAS, 1485
  * Heggie (1975) Heggie D. C., 1975, MNRAS, 173, 729
  * Heggie & Hut (2003) Heggie D., Hut P., 2003, the Gravitational Million-Body Problem: A Multidisciplinary Approach to Star Cluster Dynamics, Cambridge University Press, 2003
  * Hénon (1961) Hénon M., 1961, AnAp, 24, 369
  * Hénon (1971) Hénon M., 1971, Ap&SS, 13, 284
  * Hénon (1975) Hénon M., 1975, IAUS, 133, IAUS…69
  * Hills (1975) Hills J. G., 1975, AJ, 80, 809
  * Hurley, Pols, & Tout (2000) Hurley J. R., Pols O. R., Tout C. A., 2000, MNRAS, 315, 543
  * Hosek, et al. (2019) Hosek M. W., et al., 2019, ApJ, 870, 44
  * Lightman & Shapiro (1977) Lightman A. P., Shapiro S. L., 1977, ApJ, 211, 244
  * Lu et al. (2013) Lu J. R., Do T., Ghez A. M., Morris M. R., Yelda S., Matthews K., 2013, ApJ, 764, 155
  * Marks et al. (2012) Marks M., Kroupa P., Dabringhausen J., Pawlowski M. S., 2012, MNRAS, 422, 2246
  * Miller & Lauburg (2009) Miller M. C., Lauburg V. M., 2009, ApJ, 692, 917
  * Kroupa (2001) Kroupa P., 2001, MNRAS, 322, 231
  * Schneider, et al. (2018) Schneider F. R. N., et al., 2018, Sci, 359, 69
  * Spitzer & Hart (1971) Spitzer L., Hart M. H., 1971, ApJ, 164, 399
  * Spitzer (1987) Spitzer L., 1987, Dynamical evolution of globular clusters, Princeton University Press
  * Wang et al. (2015) Wang L., Spurzem R., Aarseth S., Nitadori K., Berczik P., Kouwenhoven M. B. N., Naab T., 2015, MNRAS, 450, 4070
  * Zhang et al. (2018) Zhang Z.-Y., Romano D., Ivison R. J., Papadopoulos P. P., Matteucci F., 2018, Natur, 558, 260