The impact of a top-heavy IMF on the formation and evolution of dark star clusters

The Spitzer instability leads to the formation of a black hole sub-system (BHSub) at the center of a star cluster providing energy to luminous stars (LSs) and increasing their rate of evaporation. When the self-depletion time of the BHSub exceeds the evaporation time of the LSs, a dark star cluster (DSC) will appear. Using the NBODY7 code, we performed a comprehensive set of direct \Nbody simulations over a wide range of initial conditions to study the pure effect of the top-heaviness of the IMF on the formation of the DSC phase. In the Galactic tidal field, top-heavy IMFs lead to the fast evaporation of LSs and the formation of DSCs. Therefore, DSCs can be present even in the outer region of the Milky Way (MW). To successfully transition to the DSC phase, the MW Globular Clusters (GCs) must possess an initial BH mass fraction of $\widetilde{\mathit{M}}_\mathrm{BH}(0)>0.05$. For star clusters with $\widetilde{\mathit{M}}_\mathrm{BH}(0)>0.08$, the DSC phase will be created for any given initial density of the cluster and Galactocentric distance. The duration of the cluster's lifetime spent in the DSC phase shows a negative (positive) correlation with the initial density, and Galactocentric distance of the star cluster if $\widetilde{\mathit{M}}_\mathrm{BH}(0)\leq 0.12$ ($\widetilde{\mathit{M}}_\mathrm{BH}(0)\geq 0.15$). Considering the canonical IMF, it is unlikely for any MW GCs to enter the DSC phase. We discuss the BH retention fraction in view of the observed properties of the GCs of the MW.


INTRODUCTION
Star clusters have been recognized as one of the key environments where black holes (BHs) are formed.However, it was initially believed that BHs experience a natal kick during their supernova explosion, which could accelerate them over the escape velocity, resulting in almost all BHs leaving the cluster as soon as they are formed.The uncertainty about the efficiency of the natal kicks has led to doubts about the number of BHs initially remaining in star clusters.Observations of Galactic low-mass X-ray binaries provide some insight into possible BH natal kicks.Integrating the trajectories of low-mass X-ray binaries containing BHs within the Milky Way (MW) reveals that while the observational parameters ★ E-mail: haghi@iasbs.ac.ir of some systems could be explained with no or small natal kicks, others are better described when considering a relatively large natal kick (Repetto et al. 2012;Repetto & Nelemans 2015;Mandel 2016;Repetto et al. 2017).
Over the last decade, there has been a significant change in our understanding of BHs in old globular clusters (GCs).Maccarone et al. (2007) identified, for the first time, a BH X-ray binary candidate inside a GC in the galaxy NGC 4472.Several more BH candidates have also been discovered in extragalactic GCs (Shih et al. 2010;Barnard et al. 2011;Maccarone & Warner 2011;Saracino et al. 2022).Strader et al. (2012) also discovered two BHs within the Galactic GC M22 by radio observations, marking the first time a MW GC had ever shown signs of having BH candidates.Considering that these BHs are accreting from white dwarf (WD) companions, and by using formation and survival rates calculated by Ivanova et al. (2010), about 5-100 BHs are thought to be present in M22.Furthermore, several BH candidates have also been discovered in other MW GCs (Chomiuk et al. 2013;Miller-Jones et al. 2015;Giesers et al. 2018).
It has been shown in several theoretical studies that GCs might actually be able to retain significant numbers of BHs (Breen & Heggie 2013;Morscher et al. 2013Morscher et al. , 2015;;Pavlík et al. 2018;Wang 2020).By investigating the evolution of two-component clusters consisting of a population of BHs co-existing within a background cluster of low-mass stars, Breen & Heggie (2013) found that the dynamical ejection rate of BHs is lower than previously thought.They used analytic calculations and direct -body simulations, which both suggest that the exchange of energy between the BHsubsystem (BHSub) and the other stars is ultimately controlled by the entire cluster.As a result, the rate of energy production in BHSub, as well as its depletion rate, is also regulated by the whole cluster.In other words, the dynamical evolution timescale of BHSub follows the evolutionary timescale of the whole cluster.They also concluded that a BHSub can survive roughly for 10 ×  rh (where  rh is the half-mass radius relaxation time).Therefore, if  rh is long enough, it is expected that MW GCs still host BHs.Breen & Heggie (2013) suggested that the collapse of the visible core occurs approximately when all BHs left the cluster.Since only about 20 per cent of the MW GCs are identified as core collapsed (Djorgovski & King 1986), it is possible that up to 80 per cent of the Galactic GCs still have populations of BHs.
Several studies, employing direct -body simulations, have shown that retaining BHs within certain clusters is necessary for reproducing their observable evidence.For instance, Mackey et al. (2008) suggested that the observed increase in core radius with age in star clusters of the Large Magellanic Cloud could be explained by a significant retention fraction of BHs because a population of heavy dark remnants inflates the core radius measured from the visible stars in the cluster (Merritt et al. 2004).Peuten et al. (2016) showed that models with a 50 per cent BH remnant fraction can reproduce the lack of mass segregation within NGC 6101.In a recent study, Gieles et al. (2021) demonstrated that the presence of a BH population within the Palomar 5 cluster, accounting for approximately 20 per cent of the cluster's present-day mass, can give rise to the extended tidal tails and large half-light radius observed in the cluster.Furthermore, Torniamenti et al. (2023) compared density profiles of direct -body models with the Gaia data of the Hyades open cluster and concluded that the observations are best reproduced by models with 2-3 BHs at present.The need to have a population of dark heavy remnants has been also put forward to interpret the central cusp in the velocity dispersion and surface brightness profiles of  Cen.Zocchi et al. (2019) fitted the velocity dispersion profile of  Cen using a modeled cluster that contained segregated BHs.Moreover, Baumgardt et al. (2019b) found that a model with 4.6 per cent of the mass of  Cen in a centrally concentrated cluster of BHs can fit all available data.
In addition to the researches focused on observational constraints of BH populations in specific clusters, recent works provide broader and more statistically robust constraints on the present-day BH content in numerous MW GCs.Various studies have been conducted using different methods, like the best-fitting multi-mass models, evaluations of visible mass segregation, and assessments of the central surface brightness of several GCs (Askar et al. 2018;Weatherford et al. 2018Weatherford et al. , 2020;;Dickson et al. 2023).These studies suggest that a moderate fraction of the present-day total mass in the form of BHs is generally necessary to explain the observations.They identified massive (> 10 5 M ⊙ ) GCs retaining especially large populations of BHs, of a total masses exceeding 10 3 M ⊙ per cluster.
As another observational example, one can refer to star cluster IRS 13E, which is an extremely compact stellar association of a few young, massive stars that are located close to the Galactic center and have survived the extreme tidal field by being bound by an invisible mass (Maillard et al. 2004) despite having relatively high-velocity (≃ 200 km s −1 ) stars (Fritz et al. 2010).Banerjee & Kroupa (2011) showed that the dark component could be an ensemble of stellar-mass BHs.They predicted that at distances close to the Galactic centre, rapid tidal stripping of star clusters by the strong tidal field can expose their BHSub, which may host a few orbiting stars.This occurs when the evaporation timescale of stars from the outer regions of the cluster is shorter than the encounterdriven self-depletion timescale of its central BHSub.Such clusters appear as highly super-virial star clusters with a large dynamical mass-to-light ( dyn /) ratio.These objects belong to a type of compact stellar populations that Banerjee & Kroupa (2011) called "dark star clusters" (hereafter DSCs).
Some observational evidence has been reported in the past few years for clusters that could be potential candidates for DSCs due to their high  dyn / ratios.For instance, Taylor et al. (2015) measured the dynamical properties of 125 compact clusters in the giant elliptical Galaxy NGC 5128.However, they could not identify any cluster with unusual kinematic properties in the dynamical mass range 10 5 <  dyn / M ⊙ < 10 6 .For  dyn / M ⊙ > 10 6 they observed two distinct sequences of star clusters in the  dyn / ∝   dyn plane.The sequence described by the slope  = 0.33 ± 0.04 exhibited  dyn / < 10 M ⊙ / ⊙ , while the sequence described by  = 0.79 ± 0.04 exhibited 10 M ⊙ / ⊙ <  dyn / < 70 M ⊙ / ⊙ , which is termed the DSC sequence.
The initial mass function (IMF) of stars within an embedded cluster is also one of the most important initial conditions that play a significant role in star cluster evolution.Most studies of resolved stellar populations in the disk of the MW showed that stars form following an IMF that has a universal form (Bastian et al. 2010;Kroupa 2001;Kroupa 2002), which is referred to as the "canonical" IMF.The shape of the stellar IMF of a star cluster near its upper mass limit is a focal topic of investigation as it determines the high-mass stellar content and hence the dynamics of the cluster in its embedded phase.Several observational and theoretical studies suggest that the IMF slope for massive stars in GCs depends on the initial cloud density and metallicity (), such that the IMF becomes increasingly top-heavy with decreasing  and increasing gas density of the forming object (Marks et al. 2012).The need for a top-heavy IMF also has been put forward to explain the observed trend of  and  dyn / ratio found in the M31 GC system, which shows a discrepancy with the stellar population synthesis (SPS) prediction (Zonoozi et al. 2016;Haghi et al. 2017).If this were to be the case, then the evolution of GCs is likely to be significantly different from the case of the canonical IMF due to the different mass-loss rates and number of BHs formed, consequently affecting our understanding of GC survival.
In this paper, we will explore the formation of DSCs for star clusters starting with a top-heavy IMF.We want to calculate a comprehensive grid of models over a wide range of initial halfmass radii ( h,i ), Galactocentric distances ( G ), , and varying the IMF-slope in the high-mass range to investigate the starting time of the DSC phase.We aim to shed light on the pure effect of the top-heaviness of the IMF as well as different values of the compact remnant retention fraction on the starting time and the lifetime of the DSC phase by means of direct -body simulations.This paper is organized as follows: In Section 2, we describe the initial setup of the -body models and our simulation method.In Section 3, we describe the dynamical process of BHSub and features of the DSC phase.The main results are presented in Section 4. Finally, Section 5 consists of a discussion and the conclusions.

DESCRIPTION OF THE MODELS
In order to quantify the effect of a top-heavy IMF on the formation and evolution of DSCs, we use the collisional -body code "NBODY7" (Aarseth 2012) which is an immediate descendant of the widely used NBODY6 direct -body evolution code (Aarseth 2003;Nitadori & Aarseth 2012).NBODY6/7 uses a fourth-order Hermite integration scheme, an individual time step algorithm to follow the orbits of cluster members, and invokes regularization schemes to deal with the internal evolution of small-N subsystems.
In addition, NBODY6 treats single and binary stellar evolution in a comprehensive way from the zero-age main sequence (ZAMS) until their remnant phases, incorporating the SSE/BSE routines and analytical fitting functions developed by Hurley et al. (2000).NBODY7 incorporates important updates on two aspects of the single stellar evolution process.Firstly, compact object masses are assigned according to the prescription from Belczynski et al. (2008).Secondly, the code implements the model presented by Vink et al. (2001) to account for mass loss due to stellar winds, as described in Belczynski et al. (2010).The high precision and efficient recipes for dealing with strong gravitational encounters, such as binary stars, are implemented in NBODY7 through the Algorithmic Regularization Chain of Mikkola & Tanikawa (1999) instead of the classic Chain Regularization in NBODY6 (Mikkola & Aarseth 1993).By utilizing this improved algorithmic regularization method for compact subsystems, it becomes possible to achieve general relativistic treatment through post-Newtonian terms and realistic parameters.Furthermore, this method provides a more thorough and reliable treatment of dynamically forming multiple systems in dense environments, particularly those involving massive objects like BHs.
A grid of 50 modeled star clusters is simulated to explore the To cover the canonical (Kroupa 2001) and top-heavy IMF (Marks et al. 2012;Kroupa et al. 2013),  3 is varied from 2.3 (canonical Salpeter value) to 1.7.We perform three sets of simulations, one with a canonical IMF (A:  3 = 2.3), and other sets with top-heavy IMFs (B:  3 = 2.0 and C:  3 = 1.7).An overview of the performed simulations is given in Table 1 and Table 2.We performed simulations with an initial cluster mass of 3 × 10 4 M ⊙ and let the models evolve for 13.2 Gyr.The dissolution time is defined to be the time when the number of stars in the cluster declines to 10 stars.The initial positions and velocities of the stars in the cluster are chosen according to a Plummer phasespace distribution function (Plummer 1911;Aarseth et al. 1974;Kroupa 2008) in virial equilibrium.Initially, the models are not mass-segregated and do not include primordial binaries.However, all types of binaries and higher multiplicity systems are allowed to form during the evolution.Some of these dynamically formed binaries are retained in the clusters over their entire evolution.
It should be noted that the BH natal kicks are, to date, poorly constrained and understood from both observational and theoretical points of view (Section 1).A common model (Belczynski et al. 2008;Fryer et al. 2012) for supernova natal kick magnitude assumes NS-like kicks (Hobbs et al. 2005) for BHs as well, but which are scaled down linearly with an increasing material fallback fraction, the so-called canonical supernova kicks, upon which mass fallback typically results in about half of BHs being ejected at birth (Chatterjee et al. 2017).However, in this paper, we assumed the extreme limit for the retained BHs in modeled clusters (i.e.full retention fraction of BHs), to examine the pure effect of the initial BH mass fraction on the evolution and formation of the DSC phase.Therefore, the natal kicks at the time of the stellar remnant formation of NSs and BHs are assumed to be negligible, allowing all of them to remain in the cluster.A detailed discussion of the retention of stellar remnants in star clusters and ultra-compact dwarf galaxies is available in Jeřábková et al. (2017) and Pavlík et al. (2018).
The initial velocity of the model star clusters is set for them to move on a circular orbit through the host Galaxy which is made up of three components: a central bulge, a disc, and a phantom dark matter halo potential (Lüghausen et al. 2015) that is scaled so the circular velocity at 8.5 kpc is 220 km/s.The bulge is modeled as a central point mass with a mass of 1.5 × 10 10 M ⊙ .The gravitational potential of the disc is represented by the Miyamoto & Nagai (1975) profile, we used values of  = 4 kpc (disc scale length) and  = 0.5 kpc (Galactic thickness), while for the disc mass, we adopted  d = 5 × 10 10 M ⊙ , as suggested by Xue et al. (2008).We adopt a logarithmic potential for the dark matter halo of the MW, here,  being the distance from the Galactic centre.The constant  c is chosen such that the combined potential of the three components yields a circular velocity of  ∞ = 220 km s −1 in the disc plane at a distance of 8.5 kpc from the Galactic centre.
All simulations were carried out on desktop workstations with Nvidia 1080 Graphics Processing Units at the Institute for Advanced Studies in Basic Sciences (IASBS).

DARK STAR CLUSTER
BHs in star clusters are formed in the first 10 Myr of the cluster's evolution with masses of about 10-100 times larger than the average stellar mass of the cluster (Belczynski et al. 2008(Belczynski et al. , 2010)).If a large number of BHs remain in the cluster, they cannot reach energy equipartition with the low-mass stars and hence undergo runaway segregation toward the centre of the cluster.The dynamical friction on the dense stellar background segregates the BHs to the cluster center (Spitzer 1987) leading to the formation of a BHSub in the central part of the cluster.This process, which is known as Spitzer's mass stratification instability, precludes thermal equilibrium.The Spitzer instability criterion is ( 2 / 1 ) ( 2 / 1 ) 3/2 > 0.16, where  1 and  2 are the average mass of light and heavy stars, respectively; and  1 and  2 denote the total mass contained in the light and heavy components, respectively.The BHs segregation towards the central part of the cluster occurs over the following timescale: where  cc is the core-collapse timescale, ⟨ BH ⟩ and ⟨⟩ are the average mass of BHs and all stars, respectively.The value of  cc depends on the initial half-mass relaxation time and is estimated to be about  cc = 15 rh for a single mass system and  cc = 0.2 rh for a realistic mass spectrum (Heggie & Hut 2003;Baumgardt et al. 2003;Fujii & Portegies Zwart 2014).The centrally segregated BHSub is dynamically very active, and many BH-BH binaries (BBHs) form via the three-body interactions in the dense stellar environment (Spitzer 1987;Heggie & Hut 2003).Indeed, BHs of different masses can interact in a single close encounter and share their kinetic energy (KE).As a result, the massive BHs in the interaction make a binary BH, while the less massive BH picks up the excess KE released in the encounter and ejects to a higher (less bound) orbit.
Through the subsequent encounters between BBHs and single BHs, the binary system becomes tighter, while the single BHs are scattered to higher orbits injecting their newly gained KE to the whole cluster via two-body interactions (Banerjee 2017), which leads to an expansion of the cluster.The scattered BHs eventually sink back to the cluster centre due to dynamical friction.With each encounter, the BBHs become more tightly bound and gain more recoil velocity.If the BBHs are sufficiently tight, during the next encounter, a significant amount of KE can be transferred to either the BBHs or single BH, leading to ejection from the cluster.This process is responsible for the self-depletion of the BH population from the clusters.
The self-depletion timescale depends on the number of BHs in the core.A cluster that retains a larger fraction of BHs will require more time for the dynamical self-emission of BHs compared to a similar cluster with fewer BHs.In other words, the dynamical selfejection of BHs is a self-regulating process in which the BH ejection times are prolonged for more massive and numerous BHs (Banerjee 2017).In star clusters orbiting in the inner part of the Galaxy, the outer skirt of star clusters strips rapidly due to the strong tidal field.When the removal timescale of low-mass stars from the outer region of the cluster due to the stronger tidal field is shorter than the selfdepletion timescale of its BHSub, a new kind of star cluster (i.e., DSCs) which is dominated by BHs can be formed.
This evolutionary phase of star clusters is predicted for the first time by Banerjee & Kroupa (2011) using numerical simulations.During this phase, the cluster includes a few low-mass stars orbiting the central BHSub that observationally appear to be super-virial with a high  dyn / ratio.This is due to the velocity dispersion of the remaining stars being enhanced by the unseen BHs.The starting time of this phase can be defined as the time at which observable luminous stars (LSs), i.e., the nuclear-burning stars and the WDs, appear to be unbound or significantly super-virial.In other words, the DSC phase starts when the virial coefficient ( = − /, where KE and PE are the kinetic and potential energy of the cluster) of the LSs,  * , becomes greater than 1 (Banerjee & Kroupa 2011).Therefore, we use the criterion of  * > 1 to determine the starting time of the DSC phase.It should be noted that all stars within twice the tidal radius are considered when calculating .
As an example, Figure 1 shows some projected snapshots of 2.0Gyr a modeled cluster (model A4) that enters into the DSC phase at about 2 Gyr.Stars with masses > 10 M ⊙ that evolve to NSs and BHs (Heger et al. 2003) are depicted in red at  = 0.In the next snapshots, as the cluster evolves, BHs (black circles) and NSs (blue circles) are formed and segregate to the centre of the cluster, while many of the low-mass stars evaporate.Due to the BH and NS segregation and frequent encounters between BBHs and BHs, the released KE injects into the LSs and causes the whole cluster to expand.Therefore, the KE of the luminous low-mass stars increases and exceeds the self-equilibrium condition ( * > 0.5).When stars in the outer region evaporate,  * rises and the luminous sub-system begins to become super-virial, where the first signs of the DSC phase appear at  = 2.0 Gyr ( * ≥ 1.0).The following snapshots show the evolution of the DSC, where the LSs and even the NSs are stripped by the external tidal field, but the BHSub survives to the final stages of cluster dissolution.Note that the spherically asymmetric extensions of the cluster seen in Figure 1, especially in the panels at  = 1.8 and 2 Gyr, are the tidal tails.The presence of the BHSub elevates the background star density by accelerating the evaporation rate of LSs.Therefore, one of the major applications of the DSC phase is the enhancement of the tidal tail/stellar stream density (Gieles et al. 2021).
Figure 2 displays the time-evolution of various dynamical parameters of model cluster A4, that can be used to introduce the DSC phase.The top panel shows a comparison of the virialcoefficient of the whole cluster's members (, blue line) enclosed by the tidal radius and the virial-coefficient of LSs ( * , red line).As expected for self-gravitational systems,  remains constant at around  = 0.5 (Heggie & Hut 2003), while  * rises above 1.0 when a significant fraction of LSs are tidally removed.The middle panel of Figure 2 illustrates that, at the onset of the DSC phase, approximately 90 per cent of LSs have escaped the cluster, while only half of the BHs have been ejected.Furthermore, it highlights that NSs evaporate earlier than BHs.Therefore, in the final stages of the cluster evolution, where the potential of the dark sub-system gradually becomes important, the cluster shows its DSC phase that is dominated by BHs.The ratio of the dynamical mass to photometric mass ( dyn / photo ) is another parameter that can indicate the DSC phase, where  photo is measured by summing the masses of all luminous objects (consisting of nuclear-burning stars and WDs).The bottom panel of Figure 2 shows that both  dyn / photo and  dyn / remain almost constant until the final stage of the evolution, increasing dramatically around the starting time of the DSC phase.
Several BBHs (and other types of binaries) form during the lifetime of our model clusters.Of these, 54 BH-BH, 3 BH-LS, 1 BH-NS, and 1 NS-LS were formed in the A4 model.Our simulations have shown that only a few BH mergers occur, all of which happened in the first 500 Myr.The number of BH mergers in GCs generally increases linearly with the cluster mass (Kremer et al. 2020).However, for GCs with top-heavy IMFs, the number of BH mergers increases super-linearly with mass (Weatherford et al. 2021).So, GCs of a typical mass should actually be significant BH merger sources, especially GCs born with top-heavy IMFs.However, the rate of BBH mergers is more intense in the early stages of cluster evolution and substantially decreases over time, especially for GCs with canonical IMFs (Banerjee et al. 2010;Abadie et al. 2011;Weatherford et al. 2021).This means that if the cluster evolves into a DSC during the later stages of its evolution, it is unlikely to detect BBH mergers that produce detectable gravitational waves during the DSC phase.Nevertheless, the remnant sub-system in DSCs can be a rich source for X-ray binaries and soft gravitational waves due to the large number of binaries with an accretor component, i.e., a BH or NS (Tauris & van den Heuvel 2006).The characteristics of the cluster that would evolve into the DSC phase can be determined as follows: • The evaporation time of LSs ought to be shorter than the BHSub self-depletion timescale.
• LSs go out of equilibrium and become super-virial.As a consequence,  * experiences a substantial increase, sometimes surpassing 100.
• Accelerating the evaporation rate of LSs enhances tidal tail/stellar stream density.
• The dynamical mass-to-light ratio of the cluster increases during the DSC phase, almost like  * .
• The ratio of BHSub mass to the total mass of the cluster, increases over time, and the value of  dyn / photo grows strongly similar to  * .At  * = 1,  dyn / photo will be approximately equal to 1.4, according to the average value of our simulated model results.
• The presence of a BHSub in a cluster can lead to the formation of a significant number of binaries containing BHs and BBHs, which can serve as sources for X-ray binaries and soft gravitational waves.If the initial cluster is very dense, BBH merging can happen, which is the source of observable gravitational waves.We expect this event to occur early in the evolution of the cluster.However, if the cluster evolves into the DSC phase at the late stages of its evolution, we would not expect to detect a gravitational wave signal in the DSC phase.
In this paper, we define the scaled DSC lifetime as follows: where  DSC is the time interval during which the cluster is in its DSC phase and  cluster is the total lifetime of the cluster.The value of  DSC determines how much of the cluster's lifetime is in the DSC phase.A zero value of  DSC implies that the cluster will not reach the DSC phase at all.On the other hand, the closer the value of  DSC is to 1, the more the cluster spends its lifetime in the DSC phase.Therefore,  DSC is a good parameter to determine how much the BHSub dominates the cluster evolution.To assess the strength of BHSub dominance, one can calculate the difference between the self-depletion time of the BHSub and the evaporation time of LSs, divided by the BHSub self-depletion time (( dep −  eva )/ dep ).A negative value of this term indicates that the cluster will not evolve to the DSC phase.This term is directly proportional to  DSC , which means that as the value of this term increases, the value of  DSC will also increase.

RESULTS
In order to study the influence of a top-heavy IMF on the formation and evolution of the DSC phase in detail, we calculated two sets of models, one in which the IMF is canonical ( 3 = 2.3), and one set with an initially top-heavy IMF ( 3 = 2.0 and 1.7).Both will be discussed in the following.The IMF of models in the following section (Section 4.1) are canonical.Then we evaluate the effect of the top-heaviness of the IMF on the results (Section 4.2).In Section 4.3, we combined the results from the previous two sections and identified the transition boundary to the DSC phase for various levels of top-heaviness.

Canonical IMF
Before we start looking into the effect of top-heavy IMF, we examine how changing  h,i ,  and  G of a star cluster with a canonical IMF influences the formation and evolution of the DSC phase.This will allow us to make inferences about the sensitivity of the results on choosing these crucial initial parameters.

The impact of the initial half-mass radii on the formation of the DSC phase
Using -body models of clusters with BHs, Gieles & Gnedin (2023) showed that the initial density is a critical parameter in setting the dynamical retention of BHs such that the mass-loss rate of the BH population due to dynamical ejections before the cluster fills the tidal radius only depends on the square root of the initial density.
In this study, we aim to investigate how the initial density of star clusters affects the formation and evolution of DSCs.Since all the modeled clusters have the same initial mass of 3 × 10 4 M ⊙ , varying the value of  h,i will result in different initial densities within the half-mass radii,  h,i = 3 cluster /8 3 h,i .We calculate models with different  h,i of 1, 3, and 5 pc, moving on a circular orbit at  G = 3 kpc.The evolution of the virial coefficient ( * ) of these clusters is shown in the top panel of Figure 3.The modeled cluster with  h,i = 5 pc (A11) is not gravitationally bound enough to resist the escape of background stars that have acquired the energy generated in BHSub.Therefore, LSs evaporate rapidly and the cluster reaches the DSC phase at about  = 1 Gyr.The model with  h,i = 3 pc (A4) reaches the DSC phase later at about  = 2 Gyr and the model with  h,i = 1 pc (A2) does not experience the DSC phase before dissolution.The bottom panel of Figure 3 shows the time evolution of the number of LSs (lime line), NSs (cyan line), and BHs (black line) for these clusters.As can be seen, the number of LSs decreases significantly in model A11 ( h,i = 5 pc), while the number of BHs remains almost constant until the final stage of the cluster evolution, when the cluster enters the DSC phase and the cluster becomes dominated by BHs.After the evaporation of LSs and NSs from the cluster, the BHSub starts to be tidally stripped.
The models with the smaller half-mass radii ( h,i = 3 pc and 1 pc pc) are relatively more gravitationally bound, which leads to slower evaporation of LSs, while the rate of few-body encounter is higher within the BHSub, leading to a faster BH self-depletion.In other words, the energy generated from the BHSub is not sufficient to bring the LSs to the escape velocity.On the other hand, collisions in the BHSub are more frequent in clusters with higher densities, leading to faster BH self-depletion compared to the LS evaporation.As a result, the number of BHs decreases with a steeper slope compared to the cluster with  h,i = 5 pc (Figure 3, bottom panel).In the case of  h,i = 1 pc, the cluster does not exhibit any sign of the DSC phase and LSs never become super-virial as the virial coefficient does not exceed  * = 1 (top panel of Figure 3, red curve).As shown in the bottom panel of Figure 3, while the number of BHs decreases with a relatively sharp slope, NSs start to replace them by accumulating in the central part as a result of the Spitzer instability.
In a nutshell, when BHs segregate due to dynamical friction, they give their initial KE to low-mass stars.In addition, during the super-elastic encounters between BBHs and single BHs, a certain amount of KE is imparted to background stars (Section 3).This amount of energy released from the BHSub can take LSs to escape velocity, thus increasing their evaporation rate, and the cluster eventually reaches the super-virial phase.In dense star clusters, the gravitational potential is deeper; hence, the released energy from the BHSub would not be enough to bring the LSs to escape velocity.So, the influence of the BHSub on the evaporation rate of LSs is more pronounced in low-density clusters compared to dense star clusters.On the other hand, the BHs are immersed in a deep central potential, before the evaporation of the LSs.During this period, the escape rate of BHs depends mostly on the few-body encounters between BBHs and single BHs rather than the tidal field of the host galaxy (Breen & Heggie 2013;Wang 2020).The LSs in the cluster halo shield the BHSub from the Galactic tidal force.Since few-body encounters in the BHSub are more frequent in a denser cluster, the BHSub will be depleted faster.As a consequence, in a dense cluster, the self-depletion time of the BHSub is shorter than the evaporation time of LSs and the cluster will not evolve into the DSC phase.
It is important to note some details about the relationship between the evaporation timescale of LSs and the initial density of clusters.When we try to determine whether the evaporation time of LSs increases or decreases with initial density, we need to consider two conflicting effects.Firstly, reducing the cluster density makes it easier for LSs to evaporate by lowering the cluster escape velocity.This can also be thought of as the cluster overflowing its tidal boundary.Secondly, lower density also slows down the evaporation of LSs by increasing the two-body relaxation timescale (Spitzer 1987;Heggie & Hut 2003).Which of these competing effects dominates depends on how tidally filling the cluster is at birth.However, the DSC phase is determined by the discrepancy between the timescales of LS evaporation and BHSub self-depletion, not by their individual rates.Indeed, if the opposite scaling is achieved (faster evaporation of LSs with higher density), dense clusters persistently fail to transition to the DSC phase.In such cases, while the evaporation time of LSs of the denser cluster becomes shorter than that of the low-density counterpart, the self-depletion time of their BHSub also contracts due to heightened encounters within BHSub, offsetting the faster evaporation time of LSs.The stronger gravitational binding of dense clusters restricts the energy injected by the BHSub to enhance the evaporation rate of their LSs.Therefore, in dense clusters, the evaporation rate of LSs is mainly determined by two-body relaxation.In clusters with a lower density, although the evaporation time is extended, the selfdepletion time is also prolonged.Since these clusters have a lower gravitational binding, the BHSub can exert a dominant influence on the evaporation rate of LSs, which facilitates the transition to the DSC phase.
The evaporation and the BHSub self-depletion timescales are compared for clusters with different densities located at  G = 3 kpc in Figure 4.The cluster evaporation time and the BHSub depletion time are defined as the times when 85 per cent of the initial mass of LSs and BHs are ejected from the cluster, respectively.The DSC phase appears when the evaporation time (red line) of LSs is shorter than the BHSub self-depletion timescale (blue line).As the value of  h,i increases, the depletion and evaporation times become closer to each other.This leads to a decrease in  DSC , indicating that the cluster spends less time in the DSC phase.As expected, for higher density, the BHSub self-depletion timescale becomes faster than the LS evaporation timescale, meaning that the DSC phase does not appear.
For a cluster with a higher stellar concentration, we expect a  faster depletion of the BHSub.However, as shown in Figure 4, for clusters located deep in the Galactic potential ( G = 3 kpc), the selfdepletion timescale of compact models (such as A2) is longer than that of extended clusters (A4 and A11).This is due to the tidal effect of the host Galaxy, which strips the BHSub after LS evaporation.Therefore, under isolated conditions, we would expect the BHSub of compact models to dissolve faster than extended clusters.

The impact of the Galactocentric distance on the formation of the DSC phase
It is well known that by increasing the Galactocentric distance of an orbiting cluster, which means reducing the tidal field of the host galaxy, the evaporation timescale increases.Conversely, the escape rate of BHs from the BHSub (i.e., the self-depletion timescale of the BHSub) is predominantly controlled by the frequency of few-body encounters within the BHSub, rather than the tidal field, as long as the halo of low-mass stars has not evaporated.Consequently, for a tidally underfilling cluster, increasing  G does not have a major effect on the self-depletion timescale of the BHSub (Breen & Heggie 2013;Wang 2020).Therefore, we expect that  G plays an important role in the formation of the DSC phase.To investigate the effect of  G on the formation of the DSC phase, we calculated models orbiting at different Galactocentric distances,  G = 2, 3, 4, 6 and 8 kpc (models A3-A7), each starting with an initial half-mass radii of  h,i = 3 pc.
As can be seen in Figure 5, the evaporation time increases with  G and becomes longer than the self-depletion time beyond 6 kpc.Then the BHSub self-depletion time becomes independent of the tidal field and does not vary as  G increases (it will eventually converge to about 8 Gyr).The bottom panel of Figure 5 shows the lifetime of the DSC phase ( DSC ) vs.  G .The value of  DSC rises up to 600 Myr at around 6 kpc.This is because after the evaporation of LSs, in a weak external field, it takes a longer time for a BHSub to deplete.As  G increases to above 6 kpc, the  DSC decreases and eventually reaches zero, which means that the cluster does not enter the DSC phase.This is because the evaporation time increases with  G , while the self-depletion time remains almost constant.The  DSC (red line) in the bottom panel of Figure 5 describes the fraction of cluster lifetimes spent in the DSC phase. DSC decreases with increasing  G such that the cluster spends 14 per cent of its life in the DSC phase at  G = 2 kpc, while it is 1 per cent at  G = 8 kpc.

The impact of the initial metallicity of the cluster on the formation of the DSC phase
Metallicity in star clusters has a significant impact on their dynamic evolution.This is because it affects the evolution of massive stars, the rate at which they lose mass through stellar winds, and the final mass of the remnants (Trani et al. 2014).Studies have shown that the retention mass fraction of BHs at birth is higher in metal- poor clusters than in metal-rich clusters (Shanahan & Gieles 2015).Moreover, in metal-poor clusters, much more massive BHs can form compared to metal-rich clusters (Vink et al. 2001;Vink & de Koter 2005;Belczynski et al. 2010).As a result, the average mass of BHs and their retained number increases as the metallicity of the cluster decreases.The energy generated from the segregation of BHs and their dynamical interactions increases with the increasing number of BHs and their average mass.Consequently, in metal-poor clusters, the energy injection from the BHSub is larger, leading to faster expansion and a higher evaporation rate of LSs.In addition, the larger number of BHs leads to a longer self-depletion time for the BHSub.Both of these factors lead metal-poor clusters to spend a larger fraction of their lifetime in the DSC phase.
To investigate the effect of  on the DSC phase, we consider three clusters with the same mass,  h,i , and  G , but different metallicities of 0.05 ⊙ , 0.25 ⊙ and  ⊙ (models A4, A8, and A9).The BHSub masses of these models were 1724 ⊙ , 1226 ⊙ , and 734 ⊙ , respectively, at birth.The time evolution of the half-mass radii,  h , of these clusters is shown in Figure 6.During the first ≈ 200 Myr, the cluster with Solar metallicity expands faster than the sub-Solar metallicity clusters due to the dominant role of stellar evolution.Metal-rich clusters lose more mass through stellar winds and supernova expulsion during the early stages (Vink et al. 2001;Vink & de Koter 2005;Schulman et al. 2012;Mapelli & Bressan 2013).However, after ≈ 200 Myr, this trend reverses, and the subsolar clusters expand more rapidly due to more energy generation from segregated BHs.As a result, metal-poor clusters expand to a larger size and dissolve earlier than metal-rich clusters (consistent with Chattopadhyay et al. 2022).The  cluster of the modelled cluster with  = 0.05 ⊙ is about half of the modelled cluster with  =  ⊙ , while  DSC is twice as large for  = 0.05 ⊙ .We also compared the mass fraction of BHs retained in these clusters when 75 per cent of their initial mass had been lost.Our findings indicate that when the metallicity is decreased from  =  ⊙ to  = 0.05 ⊙ , the mass fraction of BHs rises from 0.055 to 0.115 at a point when only 25 per cent of the initial mass is left in the clusters.Summarising, these results thus show that the DSC phase is shorter at low metalicity but comprises a longer fraction of the cluster's lifetime than at high metalicity.

The Combined effect of initial density and Galactocentric distance on the formation of the DSC phase
According to Table 1,  DSC decreases for clusters constructed based on the canonical IMF by increasing  h,i and  G .We derived the best-fitting function for  DSC in dependence of log 10 ( h,i ) and log 10 ( G ) of the modelled clusters in the following mathematical form: where  and  are the best-fitting parameters which are themselves a function of log 10 ( h,i ) as follow: Figure 7 illustrates  DSC as a function of log 10 ( h,i ) and  G . DSC increases as  G and  h,i decreases and reaches its maximum of 44 per cent for the lowest values that we adopted for  G and  h,i .Star clusters with an initial density log 10 ( h,i ) > 3.7 do not experience the DSC phase.The DSC phase within the Hubble time can only be observed in clusters with  G < 16 kpc.
According to the Marks & Kroupa (2012) relation, a cluster with a birth mass of 3 × 10 4 M ⊙ has an initial half-mass radii of 0.38 pc.It's worth noting that our simulations do not account for gas expulsion, so we have to consider the starting conditions of our models as post-gas expulsion.Baumgardt & Kroupa (2007) found that assuming a star formation efficiency of about 0.3, massive clusters expand by a factor of ≃ 3 during the gas expulsion phase.This expansion factor is almost independent of the gas removal rate.By applying this expansion factor, the post-gas expulsion density of a cluster with a mass of 3 × 10 4 M ⊙ will be about log 10 ( h,i ) = 3.4.We adopt this value as a lower limit for the initial density of the MW GCs (vertical white dashed line in Figure 7).Note that the initial density of the MW GCs is on average higher than what our model suggests.This is because the MW GCs tend to be more massive than the clusters we've used in our calculations.
Based on our analysis, it is unlikely for any MW GCs to enter the DSC phase if we consider the canonical IMF.This is because, in such clusters, the BHSub self-depletion happens at a faster rate than the evaporation of LSs.Over time, the ratio of BHSub mass to the total mass of the cluster ( M BH (), eq. 3) gradually decreases until it ultimately reaches zero.The white dashed-dotted curve separates the initial conditions of the star clusters into two groups: those that enter the DSC phase before the Hubble time (below the curve) and those that enter after (above the curve).This means that even if a cluster in the upper region could eventually evolve into the DSC phase, it would do so after the Hubble time.
The MW GCs have on average lost nearly 75 per cent of their initial stellar masses through the stellar evolution or longtime dynamical evolution (Webb & Leigh 2015;Baumgardt & Sollima 2017;Baumgardt et al. 2019a).Therefore, the BH mass fraction of the modeled clusters is plotted in Figure 8, which shows the fraction of BHs when 75 per cent of the initial mass of the clusters has been lost.For comparison with the MW GCs, the mean Galactocentric distance (Baumgardt & Vasiliev 2021;Vasiliev & Baumgardt 2021) and post-gas expulsion density of the 166 GCs are shown with white circles.To calculate the post-gas expulsion density of MW GCs, we consider 4 times their present-day dynamical mass (Baumgardt & Hilker 2018, additionally, we make use of the GC database developed by Baumgardt and Sollima1 ) as their initial stellar mass.Using the Marks & Kroupa (2012) relation, we obtain the initial half-mass radii, which is multiplied by a factor of 3 to obtain the radius of star clusters after gas expulsion (Baumgardt & Kroupa 2007).The black dashed curve corresponds to the mass fraction of BHs predicted by stellar evolution (birth fraction), which is approximately M BH (0) = 0.04.For clusters below this curve, the BH fraction increases over time, while for dense clusters located at large  G (above this curve), the BH fraction decreases over time until it reaches zero.This black dashed curve represents the boundary beyond which the cluster's transition to the DSC phase is not possible.
The comparison with MW GCs in Figure 8 shows that almost all of them are located above the dashed curve, indicating that they cannot evolve into the DSC phase (as shown in Figure 7).Therefore, assuming a canonical IMF and full initial retention of BHs, most MW GCs are nearly depleted of BHs, with only 0-1 per cent of their total mass comprising BHs. Figure 8 shows that only about 13 per cent of MW GCs could currently contain 1-4 per cent of their mass in BHs.However, it indicates that approximately 85 per cent of their BHs have escaped since their formation.These GCs are located in the inner part of the MW (mean Galactocentric distance < 4 kpc) with an initial density log 10 ( h,i ) < 4.1.
The magnitude of BH natal kick is a mystery that can be unlocked by examining the inferred present-day BH mass fraction within MW GCs.Recent studies have used different approaches to estimate this fraction, including best-fitting multimass models, evaluations of visible mass segregation, and assessment of the central surface brightness of several GCs.These studies indicate that typically, BH mass fractions ranging from 0-1 per cent of the GCs' total mass are needed to explain the observations, except for  Cen (Askar et al. 2018;Weatherford et al. 2020;Dickson et al. 2023).According to Figure 8, even if we assume zero natal kicks for BHs, only 0-1 per cent of the present-day mass of the majority of

25.
30.MW GCs (87 per cent) consist of BHs.For the rest of the MW GCs, around 85 per cent of BHs have been ejected from the cluster up to the present-day.These results highlight that achieving the presentday BH mass fraction doesn't require BHs to receive a high natal kick.Instead, even with high initial retention of BHs, a substantial number of them are depleted through few-body encounters in the core of the GCs, shaping the present-day BH mass fraction.According to Gieles & Gnedin (2023), the reduction of M BH () due to dynamical ejections depends only on the square root of the initial density relative to the tidal density before the cluster fills the tidal radius, and not on the initial relaxation timescale.Therefore, although increasing the initial modeled cluster mass primarily influences parameters such as the cluster relaxation time, it alone is unlikely to have a significant effect on M BH (), as the M BH () of the cluster is mostly determined by cluster density.As a result, even though the data in Figure 8 is specific to the modeled clusters with a mass of 3 × 10 4 M ⊙ , it could be extrapolated to more massive clusters similar to those of the MW GCs.
Note that Figures 7 and 8 demonstrate the  DSC and fraction of BHs for clusters that initially retain all of their BHs.However, if we apply natal velocity kicks, such as canonical supernova kicks with mass fallback to the modeled clusters, leads to a re-scaling of the colors in these figures.Therefore, it is more accurate to state that Figures 7 and 8 essentially depict  DSC and the fraction of BHs for clusters with M BH (0) = 0.04, irrespective of their IMF and the natal kick received by their BHs.
According to our simulations of clusters with a canonical IMF, when the cluster evolves into the DSC phase,  dyn / photo is approximately equal to 1.25, and its evolution is similar to  * and rises sharply in the DSC phase.Similarly, the dynamical mass-tolight ratio has the same evolution as  * .

Top-heavy IMF
The number of BHs formed in a cluster is determined by the highmass end of the stellar IMF.If the IMF is top-heavy, relatively more massive stars will be formed during star formation, resulting in a relatively greater number of BHs.This increase in BHs causes BH heating through mass segregation and consequently, a higher rate of few-body encounters in the BHSub, leading to a faster supervirial phase and increased evaporation rate of LSs.As a result, the dissolution timescale decreases with an increase in the slope of the mass function in high-mass stars, which is measured by the value of  3 (Chernoff & Weinberg 1990;Chatterjee et al. 2017;Haghi et al. 2020;Weatherford et al. 2021).To examine the effects of the IMF on the evolution of the DSC phase, we set up a series of models (sets B and C) with varying degrees of top-heaviness.
Figure 9 depicts the time-evolution of the number of LSs (lime line), NSs (cyan line), and BHs (black line) in clusters with different values of  3 that orbit circularly at a radius of  G = 8 kpc.For a cluster with a canonical IMF (right panel), the mass fraction of BHs at birth is M BH (0) = 0.04.In this cluster, the energy produced by the BHSub is not significant enough to affect the evaporation rate of LSs, so the BHSub dissolves before the LSs have enough time to evaporate.As the BHSub self-depletes, the NSs form a segregated subsystem located at the centre of the cluster.As the cluster approaches complete dissolution, a short-lived DSC phase appears with  DSC = 0.016.
For a cluster with  3 = 2.0 (middle panel), the mass fraction of BHs at birth is M BH (0) = 0.09.In this cluster, the BHSub produces enough energy to significantly increase the evaporation rate of LSs.Therefore, BHs remain in the cluster until the final stage of cluster evolution, while LSs and NSs have already evaporated.Consequently, the cluster spends approximately half of its lifetime in the DSC phase (  DSC = 0.555).Assuming a very top-heavy IMF with  3 = 1.7 (left panel), the cluster becomes highly dominated by BHs ( M BH (0) = 0.15), and the produced BHSub remains largely intact until dissolution, while other components have already escaped.In this cluster, which contains a significant number of BHs, the BHSub produces so much energy that LSs become super-virial already after approximately 100 Myr.This cluster is in the DSC phase for almost its entire life (  DSC = 0.958).
As illustrated in Figure 10, for the examined top-heavy models,  3 = 2.0 (B1-B9) and  3 = 1.7 (C5-C14), the self-depletion time is longer than the evaporation time at all Galactocentric distances.This means that all computed models can evolve to the DSC phase.For the  3 = 2.0 models (Figure 10 top panel), the difference between self-depletion time and evaporation time increases as  G increases, while the parameters ( dep − eva )/ dep and  DSC decrease slowly as  G increases.It's worth noting that for clusters located above ≈ 50 kpc, the self-depletion time becomes longer than the Hubble time (shown as a dashed gray line in Figure 10).Clusters evolving at orbital radii larger than 80 kpc experience neither the LSs evaporation nor the DSC phase before the Hubble time.
The modeled clusters C5-C14 with  3 = 1.7 generate intense energy from BH segregation, causing the LSs to become supervirial within approximately 100 Myr and entering the DSC phase at the same time as the BHs segregate.The tidal field does not have a significant effect at this time.Even an isolated cluster (C16) enters the DSC phase around the same time.As the LSs become super-virial, they quickly evaporate from the cluster, and increasing  G has little effect on reducing the evaporation rate.The bottom panel of Figure 10   LSs escape, only the BHSub remains, with a longer dissolution time in weaker tidal fields (larger  G ). Therefore both parameters, ( dep −  eva )/ dep and  DSC , increase with increasing  G .
Analyzing the calculated models with various densities and Galactocentric distances, we discovered the most accurate linear relationship between  DSC and log 10 ( h,i ) as well as log 10 ( G ) (as seen in Equation 7).This relationship is depicted in Figure 11 for both  3 = 2.0 (top panel) and  3 = 1.7 (bottom panel).From The results indicate that, when a cluster transitions into the DSC phase ( * = 1), the ratio of  dyn / photo is roughly 1.35 and 1.5 for  3 = 2.0 and  3 = 1.7, respectively.Clusters with a top-heavy IMF exhibit a softer evolution of  * ,  dyn / photo , and  dyn / during the DSC phase, unlike the sharp increase observed for the canonical IMF.

The DSC transition boundary for different degrees of top-heaviness
Based on the findings presented in Section 4.1 and Section 4.2, it is evident that the value of  DSC relies on several factors, including ,  h,i ,  G and  3 (i.e., M BH (0)).The larger the value of  DSC , the longer the cluster will remain in the BH-dominated state, leading to a more rapid increase in the M BH () ratio, which can eventually reach 1.0.Hence, if  DSC > 0, the M BH () ratio will continue to increase over time, and the transition boundary is defined by  DSC = 0, below which value M BH () decreases over time, ultimately resulting in the complete depletion of BHs in the cluster before the cluster dissolves.
The transition boundary at which  DSC is zero for a given M BH (0) in  G −  h,i space is plotted in Figure 12.The contours show the necessary boundary conditions for a cluster to evolve into the DSC phase, based on a specific M BH (0) value.When the value of M BH (0) increases, the region in the  G and log 10 ( h,i ) space where clusters can transform into the DSC phase expands.Figure 12 shows that for MW GCs to evolve into the DSC phase, the retention fraction of BHs at birth must be greater than 0.05.Additionally, it can be concluded that the DSC phase will occur in all values of  h,i and  G when M BH (0) > 0.08.It is important to note that the limit of the Hubble time was not accounted for in Figure 12.The determination of transition boundaries depends only on the initial mass fraction of BHs in the cluster, regardless of the magnitude of the natal kick received by them.It is expected that these transition boundaries will still be applicable for clusters that retain the mass fraction of BHs, M BH (0), even after experiencing the natal velocity kick of BHs.It is important to note that the region of space leading to a DSC is highly sensitive to M BH (0).Even canonical supernova kicks with mass fallback, let alone higher-speed kicks without fallback, reduce M BH (0) by about a factor of 2, which can alter the interpretation of the plot.This means that even realistic MW GCs with slightly top-heavy IMFs may not be capable of evolving to a DSC phase.Recently, Zocchi et al. (2019) fitted dynamical models to  Cen data and showed that models with 5 per cent of their mass in  Wang (2020) showed that in tidally filling models, M BH () increases rapidly in models with M BH (0) > 0.07, leading to the formation of BH-dominated dark clusters.However, the value of 0.07 is not constant and depends on the ratio of the average mass of BHs to the average mass of other stars and the tidal radius divided by the half-mass radius.In our study, we performed 50 simulations that accounted for stellar evolution and varied initial conditions such as M BH (0),  h,i , and  G .Our results indicate that the transition boundary for a cluster to enter the DSC phase is approximately M BH (0) = 0.08.
According to the results obtained in the previous sections, the dependency of  DSC on  G and  h,i as a function of M BH (0) can be obtained.The solid line in Figure 13  in Figure 13 separates the region in which   DSC / G is positive (II and III) or negative (I).In between (II) is the region in which   DSC / G > 0 and   DSC /  h,i < 0. As can be seen, for star clusters with M BH (0) ≤ 0.12,  DSC decreases with increasing  G and  h,i , while it increases for star clusters with M BH (0) ≥ 0.15.For clusters in this region, the energy generated from the segregation of BHs is so intense that the LSs become super-virial within nearly the first 100 Myr.Therefore, clusters with M BH (0) ≥ 0.15 spend their entire lifetime in the DSC phase.

CONCLUSION
Using numerical simulations of star clusters with a canonical IMF, Banerjee & Kroupa (2011) predicted, for the first time, the formation of dark star clusters (DSC) as a result of BH segregation to the centre of the cluster and rapid removal of stars from the outer parts of a cluster caused by the strong tidal field of the host galaxy.As a generalization of that work, we explore the formation of DSCs in star clusters starting with a top-heavy IMF.In this work, we have carried out a series of direct -body simulations of star clusters over a wide range of half-mass radii, metallicities, Galactocentric distances and IMF slope in the high-mass range to investigate the starting time and duration of the DSC phase (see Table 1 and Table 2).In all simulations, we assumed zero natal kicks for NSs and BHs.With the complete retention fraction of the BHs, the cluster becomes Spitzer unstable, which leads to the formation of a black hole subsystem (BHSub) in the central part of the cluster.
The energy produced by the BHSub causes the LSs to enter a super-virial phase and speeds up their evaporation rate.Whether or not the cluster can reach the DSC phase depends on the time between the depletion of the BHSub and the evaporation of the LSs.We use the scaled DSC lifetime (  DSC ), which is defined as the duration of the cluster's DSC phase divided by its overall lifetime, to measure how long the BHSub dominates the evolution of the cluster.We examined the dependency of  DSC on ,  h,i ,  G and  3 (or equivalently the mass fraction of the BHs at birth, M BH (0)).The main outcomes of our study can be summarized as follows: • The ratio of the dynamical mass (obtained from velocity dispersion) to the photometric mass ( dyn / photo ) and the dynamical mass-to-light ratio ( dyn /) of the cluster are other parameters that can indicate the DSC phase.These observable parameters show a strong similarity to the  * parameter and undergo a similar change.We showed that at the time when the cluster evolves into the DSC phase (i.e.,  * = 1),  dyn / photo is approximately equal to 1.25, 1.35, and 1.5 for  3 = 2.3, 2.0 and 1.7, respectively.Generally, a cluster is in the DSC phase if its dark remnant mass fraction is greater than 28 per cent assuming the WDs are luminous masses.
• If the scaled DSC lifetime is equal to zero, then the ratio of the mass of the BHSub to the total mass of the cluster, M BH (), decreases over time.This will result in the BHSub being depleted in the cluster.However, if  DSC > 0, then the M BH () ratio increases and eventually reaches 1.0.
• The presence of a BHSub leads to the formation of a considerable amount of binaries containing BHs and BBHs.These are excellent sources of X-ray binaries and soft gravitational waves.The merging of BBHs, which results in detectable gravitational waves, mainly occurs during the initial stages of cluster evolution.However, if the cluster evolves into a DSC during the later stages of its evolution, it is unlikely to detect such a gravitational wave event during the DSC phase.
• The lifetimes and scaled DSC lifetimes of the metal-rich ( =  ⊙ ) clusters are approximately two and a half times longer than those of the metal-poor ( = 0.05 ⊙ ) clusters.
• In order for the MW GCs to evolve into the DSC phase, the retention fraction of the BHs at birth must be greater than 0.05.We have determined the minimum value of the initial BH mass fraction, M BH (0), that guarantees the cluster will reach the DSC phase.If this fraction is greater than 0.08, the DSC phase will occur at every Galactocentric distance and with varying initial density.
• If MW GCs have followed the canonical IMF and even if their BHs have not received any natal kicks, most of them are nearly depleted of BHs at present, with only 0-1 per cent of their total mass attributed to BHs.However, around 13 per cent of MW GCs could still have 1-4 per cent of their mass in BHs.Nevertheless, approximately 85 per cent of their BHs have escaped since their birth.These GCs are located near the bulge, within a mean Galactocentric distance of less than 4 kpc, and have a low initial density of less than log 10 ( h,i ) < 4.1.Recent studies have shown that BH mass fractions ranging from 0-1 per cent of the total masses of MW GCs are typically required to explain the observations.Based on this, we can conclude that achieving the present-day BH mass fraction does not require BHs to receive a high natal kick.Even with high initial retention of BHs, a substantial number of them are depleted through few-body encounters in the core of the GCs, shaping the present-day BH mass fraction.
• The identification of a BH mass fraction exceeding 1 per cent within MW GCs that orbit the Galactic center at mean Galactocentric distances larger than 4 kpc, or those characterized by high initial densities (log 10 ( h,i ) > 4.1), can be interpreted as evidence for them having been born with a top-heavy IMF.The IMF of a GC can thus be constrained by determining the mass fraction of its BHs.
• Several studies indicate that a model comprising over 5 per cent of the mass of  Cen in a cluster of BHs that are concentrated in the center could replicate the velocity dispersion profile of  Cen without requiring an IMBH at the center of the cluster.We showed that this suggests that the IMF of  Cen may have been top-heavy in order to have more than approximately 5 per cent of its mass in a BHSub after 12 Gyr.
• We showed that the scaled DSC lifetime,  DSC , decreases with increasing the Galactocentric distance and initial density of star clusters if M BH (0) ≤ 0.12, while  DSC shows a positive correlation with both  G and  h,i if M BH (0) ≥ 0.15.

Figure 1 .
Figure 1.Top-down view of model A4.The snapshots are depicted in a frame aligned with Galactocentric coordinates.Upper left: the initial distribution of low-mass stars (dots) and massive stars (> 10 M ⊙ ) that end up as NSs or BHs (circles).Upper centre: the distribution of NSs (blue circles) and BHs (black circles) at  = 100 Myr.Upper right: the cluster starts to exceed the virial equilibrium.Lower left: the beginning of the DSC phase.Lower centre: the cluster is significantly super-virial with  * = 5.0.Lower right: the final stage of the cluster's evolution, where the number of BHs is comparable to the luminous star population.

Figure 3 .
Figure 3. Top: time variation of  * for clusters with three different half-mass radii of  h,i = 1, 3, and 5 pc (red, orange, and green lines, respectively).Bottom: the number of stars (lime line), neutron stars (cyan line), and black holes (black line) as a function of time for clusters depicted in the top panel.

Figure 4 .
Figure 4.For simulations in set A with  G = 3 kpc, the blue solid and red dashed lines represent the self-depletion and evaporation times, respectively, of BHSub and LSs.

Figure 5 .
Figure 5. Top: evaporation time of LSs (red diamonds) and self-depletion time of the BHSub (blue squares) at different  G for models A3 to A7. Fitted curves, represented by a blue solid line and a red dashed line, respectively, capture trends in self-depletion and evaporation time.The horizontal dashed line shows the Hubble time (13.2Gyr).Bottom: the DSC lifetime (the left axis) and the scaled DSC lifetime (  DSC , Eq. 6, the right axis) for different  G are shown by blue solid and red dashed lines, respectively.

Figure 6 .
Figure 6.The evolution of half-mass radii for three clusters with three different metallicities of  = 0.001 (blue line),  = 0.005 (red line), and  = 0.02 (green line).The time intervals ( DSC ) during which the clusters are in their DSC phase are shown.

Figure 7 .
Figure 7.  DSC is shown as a function of log 10 ( h,i ) and  G .The right side of the perpendicular with a dashed line shows an acceptable range for the initial density of the MW GCs.The dashed-dotted line represents the curve where the starting time of DSC equals the Hubble time.

Figure 8 .
Figure 8. BH mass fraction of the modeled clusters with a canonical IMF at the moment when they have lost 75 per cent of their mass since formation.The dashed contour represents a BH mass fraction of 0.04.The white circles indicate where the 166 MW GCs are located in the 2D space of initial density and Galactocentric distance.The vertical with dashed line is as in Fig. 7.

Figure 9 .
Figure 9.The variation of the luminous star, neutron star, and black hole populations in clusters with different slopes of  3 = 1.7, 2.0, and 2.3.

Figure 11 ,
Figure11, we can conclude that with a top-heavy IMF, many MW GCs can enter the DSC phase (for  3 = 2.0, only clusters with  G < 35 kpc) before the Hubble time.A cluster's relaxation time decreases with a reduced  h,i , leading to a decrease in segregation time (Equation4).Reduction in the segregation time of BHs can hasten the attainment of the supervirial state of LSs for clusters with  3 = 1.7.For example, the DSC phase starts for C3, C11, and C19 ( G = 16 kpc) at 53, 105, and 174 Myr, respectively.As the clusters reach the DSC phase more quickly, the  DSC increases.As a result, for clusters with  3 = 1.7, the  DSC exhibits a positive correlation with  G and  h,i , unlike as for the canonical IMF and top-heavy IMF with  3 = 2.0.The bottom panel of Figure11highlights that the clusters with an initial density of MW GCs spend almost their entire lifetime in the DSC phase.The results indicate that, when a cluster transitions into the DSC phase ( * = 1), the ratio of  dyn / photo is roughly 1.35 and 1.5 for  3 = 2.0 and  3 = 1.7, respectively.Clusters with a top-heavy IMF exhibit a softer evolution of  * ,  dyn / photo , and  dyn / during the DSC phase, unlike the sharp increase observed for the canonical IMF.

Figure 12 .
Figure 12.The contours where the  DSC is zero at M BH (0) = 0.04, 0.05, 0.06, 0.07, and 0.08 in the 2D space of  h,i and  G .Star clusters located below the curves can evolve to the DSC phase.

]Figure 13 .
Figure 13.To determine the positive or negative correlation of  DSC , the phase space is divided into three regions based on the variation of  G ,  h,i , and M BH (0) =  BH (0)/ (0).The black curve indicates the points where   DSC / h,i becomes zero for specific values of  G and M BH (0).The dashed curve represents the points where   DSC / G becomes zero for specific values of  h,i and M BH (0).

Table 1 .
Parameters of the various  -body models with different initial half-mass radii ranging from  h,i =1 to 5 pc, different Galactocentric radii ranging from G =2 to 16 kpc, and metallicity of  = (0.05, 0.25, 1)  ⊙ .The table also lists the name and main defining property of each model, the DSC lifetime (column 6), the cluster lifetime (column 7), and their ratio (column 8).The canonical IMF ( 3 = 2.3) is used in all models of set A.