New Parameters for Star Cluster Dynamics: The Role of Clusters’ Initial Conditions

We recently introduced three new parameters that describe the shape of the normalized cumulative radial distribution (nCRD) of the innermost stars in globular clusters (GCs) and trace the clusters’ dynamical evolution. Here, we extend our previous investigations to the case of a large set of Monte Carlo simulations of GCs, started from a broad range of initial conditions. All the models are analyzed at the same age of 13 Gyr when they have reached different evolutionary phases. The sample of models is well representative of the structural properties of the observed population of Galactic GCs. We confirm that the three nCRD parameters are powerful tools to distinguish systems in the early stages of dynamical evolution from those that have already experienced core collapse. They might also help disentangle clusters hosting a low-mass intermediate-mass black hole of a few hundred solar masses from cases with large concentrations of dark remnants in their centers. With respect to other dynamical indicators, the nCRD parameters offer the advantage of being fully empirical and easier to measure from observational data.


INTRODUCTION
Globular clusters (GCs) are among the densest stellar environments in the universe, comprising up to a few millions of stars in approximately spherical configurations with typical values of the half-mass radius of just a few parsecs.Various internal dynamical processes and the interaction with the external tidal field of the host galaxy shape the present-day dynamical properties of GCs and their stellar content.Two-body relaxation, in particular, plays a key role in driving GC evolution, and the manifestations of its effects include segregation of massive stars in the cluster's inner regions, and evolution towards energy equipartition.One of the key features of a cluster's internal dynamical evolution is the progressive contraction of the core as the system approaches the core collapse (CC) phase, which is eventually halted by the energy generated by binary stars (see, e.g., Heggie & Hut 2003).If stellar mass black holes (BHs) are retained within the cluster, the BH sub-system experiences an early contraction just 10-100 Myr after formation (Breen & Heggie 2013).This delays the CC of the observable stars until almost all BHs are ejected from the system.The cluster central density and central structure, and its post-CC evolution depend on the fraction of primordial or dynamically formed binaries, as well as the fraction of dark stellar remnants.The post-CC phase may be characterized by large fluctuations in the core properties driven by gravothermal effects (see, e.g., Makino 1996 and references therein).The structural evolution of a GC depends on the interplay among a number of factors including, for example, the fraction of primordial binary stars, the fraction of BHs retained in the cluster, the strength of the external tidal field, and the cluster's mass, size, and internal kinematics (see, e.g., Chernoff & Weinberg 1990;Mackey et al. 2007;Giersz et al. 2013;Breen et al. 2017;Askar et al. 2018;Kremer et al. 2018Kremer et al. , 2021;;Gieles et al. 2021;Pavlík & Vesperini 2022).The high cluster densities, especially during the CC and post-CC phases, may also enhance the frequency of close encounters and physical collisions between stars and binaries inside the core, making them efficient factories for the production of observable astrophysical exotica, such as blue straggler stars, cataclysmic variables, millisecond pulsars and tight binaries Bhat et al. with compact degenerate companions, including binary BHs that may lead to gravitational wave events (see, e.g., Rodriguez et al. 2015Rodriguez et al. , 2016;;Askar et al. 2017;Hong et al. 2017Hong et al. , 2018;;Ye et al. 2020;Antonini & Gieles 2020).The degree of mass segregation and energy equipartition also strongly depends on the cluster's dynamical phase, its stellar content, and the fraction of dark remnants (e.g., stellar-mass and intermediate-mass BHs; see, e.g., Gill et al. 2008;Alessandrini et al. 2016;Peuten et al. 2016;Aros & Vesperini 2023).The dynamical evolutionary stage of any observed GCs is thus influenced by the combination of various effects, and GCs with the same chronological ages typically have different dynamical ages, some having already reached CC, some other still being in early dynamical stages.
The observational identification of the dynamical phase of a GC and the link between its dynamical history and the different properties affecting its evolution can be quite challenging.Typically, the presence of a central power-law cusp in the observed density profile of GCs in comparison with the flat core behavior of King (1966) models is used as an indicator of CC or post-CC clusters (Djorgovski & King 1984;Chernoff & Djorgovski 1989).This diagnostic, however, can not fully capture and identify the various possible dynamical phases of a GC and the degree to which dynamical processes have altered the cluster's internal dynamical properties and stellar content.For this reason additional indicators have been proposed in recent years.These indicators are either based on peculiar populations of heavy stars (e.g., blue straggler stars) that are sensitive tracers of the dynamical friction efficiency (Ferraro et al. 2012(Ferraro et al. , 2018(Ferraro et al. , 2019(Ferraro et al. , 2020(Ferraro et al. , 2023;;Lanzoni et al. 2016), or based on the radial variations of the stellar mass function, kinematics, or distribution of stars in the cluster (Baumgardt & Makino 2003;Webb & Vesperini 2017;Tiongco et al. 2016;Bianchini et al. 2016;Bianchini et al. 2018;Leveque et al. 2021;Bhat et al. 2022).
In a recent study (Bhat et al. 2022, hereafter B22) we defined three new diagnostics of dynamical evolution based on the normalized cumulative radial distribution (nCRD) of "normal" cluster stars, i.e., stars observed at about the main-sequence turnoff level and above.The study has been performed by using Monte Carlo simulations of star clusters, and it has been further extended in Bhat et al. (2023, hereafter B23) to investigate the effects of different populations of primordial binary systems and stellar-mass BHs.In those studies we found that the proposed new diagnostics can efficiently distinguish clusters in various evolutionary phases, including CC and post-CC, and we established a connection between them and the various parameters (e.g., primordial binary fraction, content of stellar-mass BHs) affecting the dynamics of GCs.In those initial investigations we have followed the time evolution of the new diagnostics measured at different times for a few representative systems.In this paper we extend our investigation to explore the evolution of GCs starting from a much broader range of initial conditions.Instead of following an individual cluster at different chronological ages, here we study the new parameters introduced in B22 in a sample of coeval clusters that, due to their different initial conditions and evolutionary history, after 13 Gyr of evolution reached different dynamical phases.The dataset we will focus our attention on is thus similar to the observational population of Galactic GCs, which have comparably old chronological ages (of ∼ 13 Gyr), but different dynamical ages.
The paper is organized as follows.In Section 2, we describe the initial conditions of the Monte Carlo simulation survey and the methodology adopted in the following analysis.In Section 3, we discuss the time evolution of the 1% Lagrangian radius of the extracted snapshots, we show the comparison between the structural parameters of the simulated clusters and those observed in the Galactic GC population, we presents the nCRDs of survey simulations, and recall the definition of the three dynamical indicators.In Section 4 we describe the reference boundaries and models adopted for the analysis of the nCRD parameters, elaborate on the impact of the initial conditions, and discuss the ability of new diagnostics to distinguish star clusters on the basis of their dynamical evolutionary stage.Finally, the summary and conclusion of the work are provided in Section 4.

METHODS AND INITIAL CONDITIONS
The GCs analyzed in this work have been simulated through the MOCCA Monte Carlo code which includes the effects of two-body relaxation, a tidal truncation, binary-single and binary-binary encounters and, following the prescriptions based on the SSE and BSE codes by Hurley et al. (2000Hurley et al. ( , 2002)), also the effects of binary and single stellar evolution (Hypki & Giersz 2013;Giersz et al. 2013).Massive star winds and BH masses are obtained from the standard SSE/BSE prescriptions.The wind mass loss prescriptions on the MS are based on Nieuwenhuijzen & de Jager (1990) and Humphreys & Davidson (1994) for mass loss in luminous blue variable stars.We run a total of 108 simulations started from different combinations of initial conditions for the density profile (for which we adopted King 1966 models with different values of the central dimensionless potential W 0 ), number of stars N (defined as the sum of the number of single and binary stars, N = N s +N b ), tidal filling factor (TF, defined as the ratio of the half-mass to the tidal radius, r h /r t ), and galactocentric distance R g .In all cases, the metallicity is Z = 10 −3 , and the primordial binary fraction, defined as the ratio N b /(N s + N b ), is equal to 10%.As in B22, the binary properties are set according to the eigenevolution procedure outlined in Kroupa (1995) and Kroupa et al. (2013).The stellar masses range between 0.1 and 100M ⊙ following a Kroupa (2001) initial mass function.Supernova natal kick velocities for BHs are assigned either according to a Maxwellian distribution with velocity dispersion of 265 km s −1 (Hobbs et al. 2005), or according to the mass fallback procedure described by Belczynski et al. (2002), where BHs have a reduced kick velocity, hence a higher retention fraction (Arca Sedda et al. 2018;Askar et al. 2018).We use k = 0 and k = 1 to flag models adopting kick velocities following the Maxwellian distribution from Hobbs et al. ( 2005) and with reduced kick velocity for BHs, respectively.Remnant BH masses range from 5M ⊙ to up to 20M ⊙ (which is maximum BH mass from single star evolution at Z=0.001 in standard SSE).For models in which mass fallback prescriptions from Belczynski et al. (2002) are enabled, the maximum BH mass from single stellar evolution is about ∼30M ⊙ .After 50 Myr the fraction of BHs retained in the clusters in simulations with fallback on is about 45-50%.The initial values adopted for the different parameters are summarized in Table 1.Their combination corresponds to a total of 108 runs.To identify each of these simulations, we introduce a naming convention based on their initial conditions.For instance, N500 K0 RG2 TF0.1 W05 represents the simulation run with N=500K, k = 0, R g =2 kpc, TF=0.1 and W 0 =5.Out of 108 simulated GCs, only 92 survived after 13 Gyr of evolution.As in the previous studies (B22, B23), they have been analyzed from the point of view of an observer, meaning that standard procedures and approximations adopted in observational works have been applied: in each snapshot, the simulated cluster is projected onto a 2D plane, and a distance of 10 kpc from the observer has been assumed to transform the distances from the centre of the system from parsecs to arcseconds.In addition the magnitude of binary systems is computed by summing up the luminosities of both components in analogy to what is done in observed GCs where they can not be resolved.

Time evolution of the 1% Lagrangian radius
Figures 1 and 2 show the time evolution of the 1% Lagrangian radius (r 1% , i.e., the radius including 1% of the total cluster mass) for a few selected runs representative of different dynamical evolutionary histories found  tem undergoing CC1 before 13 Gyr.Conversely, the upper right panel illustrates the case of a system that experiences CC slightly after 13 Gyr.As discussed in Section 1, the simulations with k=1 where BHs are initially retained the clusters are characterized by an early inner contraction due to the early segregation and core collapse of the BH subsystem.The two lower panels of the same figure show the cases of simulations in which r 1% experiences an extended contraction phase and, after 13 Gyr of evolution, are still far from CC (which likely occurs significantly later than 15 Gyr).Fig. 2 shows representative cases of simulated clusters that, at an age of 13 Gyr, host an intermediatemass BH (IMBH) or a large population of stellar-mass BHs.The leftmost upper panel shows r 1% undergoing a contraction phase in the first 6 Gyr, followed by a rapid and significant expansion.As can be seen in the panel below, this behavior is due to the presence of an IMBH, which rapidly starts growing in mass around 6 Gyr, reaching a value of M∼ 229M ⊙ at 13 Gyr (see Giersz et al. 2015 for a detailed study of possible evolutionary paths leading to the formation of IMBHs).The IMBH becomes a dominant energy source in the core, sustaining the cluster against contraction and triggering the expansion of r 1% .The panels in the second column of Fig. 2 show a different case of simulations with an IMBH.Here, the mass growth of the IMBH is more gradual, starting already from the very beginning of the run (bottom panel).Correspondingly, r 1% shows an almost constant, slightly expanding, evolution in time (upper panel).
The remaining panels of the same figure present simulations with reduced BH kick velocity (k = 1) that both retain a large number of stellar-mass BHs for many Gyrs, but show distinct behaviors in r 1% .The top panel in the third column shows a run where r 1% is expanding for the first 7 Gyr and then undergoes contraction due to two-body relaxation (similar to the simulation named "DRr" in B23).As shown in the corresponding bottom panel, the system loses stellar-mass BHs at considerably high rate, with no BHs left at an age of 10 Gyr.The rightmost column of Fig. 2 shows a run with k = 1 but where r 1% experiences a continuous expansion for the entire cluster lifetime.The time evolution of the number of BHs in this system (rightmost bottom panel) indicates that the BH ejection proceeds at a much lower rate than in the previous case (third column in the Figure 2), and ∼ 400 such objects are still present at an age of 13 Gyr.The large majority of simulations started with k = 1 shows behaviors of r 1% like those plotted in the rightmost panels of this figure.In only two instances (simulations N500 K1 RG2 TF0.025 W05 and N500 K1 RG2 TF0.025 W07; the latter is the run named "DRe" in B23), the BH ejection rate is rapid enough for the systems to undergo a substantial dynamical evolution and reach CC and post-CC phases.
Altogether, the time dependence of r 1% within the survey simulations reveals a large variety of dynamical evolutionary histories.Based on the behavior of r 1% , we classify the simulated cluster extracted at t = 13 Gyr into three main categories: "near to CC (nCC)" if the cluster is reaching CC in the next 1-2 Gyr (as in the topright panel of Fig. 1), "CC" if they are in CC or post-CC phases (as in the top-left panel of Fig. 1), and "pre-CC" in all the other cases.The snapshots showing a rapid growth of an IMBH in their center at the time of CC (as the one shown in the upper-left panel of Fig. 2) are classified as pre-CC clusters because the substantial reexpansion they suffer brings them to have, at an age of 13 Gyr, the typical structure of poorly evolved systems, with large values of the 1% Lagrangian radius.

Comparison with Galactic GCs
To verify whether the synthetic clusters are representative of the population of GCs observed in our Galaxy, here we compare their structural parameters with those listed in the Harris (1996Harris ( , 2010 version) version) catalog.
To determine the structural parameters of the 92 simulated clusters extracted at 13 Gyr from our Monte Carlo survey, we follow the same procedure described in B22 and adopted in observational works (e.g., Miocchi et al. 2013;Lanzoni et al. 2019, see also B22).This consists in determining the projected density profile from resolved counts of stars brighter than 1 magnitude below the main sequence turnoff point (i.e., with V < V T O +1), in a series of concentric rings around the cluster center, and then search for the best-fit King model through a χ 2 minimization approach (see also Section 3.1 in B22 and Section 3 in B23).We find that the King model family fits quite well the density profiles of most (∼ 80%) the snapshots, with the exclusion of the ones that are close to and beyond CC, when a density cusp develops in the center.As discussed in B22 and B23, for the sake of homogeneity and reproducibility, and to avoid arbitrariness, the cusps are not excluded from the King fits.These density cusps are shallow compared to those presented in B22, consistent with the fact the present runs started with 10% binary fraction and in some cases retain a significant fraction of stellar-mass BHs (see also B23).From the best-fit King model to each simulated cluster, we obtained structural parameters as the core radius R c (i.e., the radius at which the central pro-jected density is halved), the half-light radius R hl (i.e., the radius including half the total projected light), and the concentration parameter c (which is defined as the logarithm of the ratio between the tidal and the King radii).We also determined the core relaxation time (t rc ) of each system following the same procedure adopted in observational studies, i.e., by using equation (10) in Djorgovski (1993) and adopting an average stellar mass < m * >= 0.3M ⊙ and a V −band mass-to-light ratio M/L V = 3.
Figure 3 shows the resulting trends among R hl , R c , c, and t rc , for the simulated clusters (empty symbols), compared to the values listed in the Harris catalog for the observed Galactic GCs (filled circles).Overall the simulations occupy the same regions of the parameter space covered by the observations.The only exception concerns the regions occupied by the Galactic GCs classified as CC and post-CC in the Harris compilation (gray filled circles in the figure ), where no simulated systems are found.This discrepancy is due to the different approach adopted for the determination of the structural parameters of these clusters, which show a central cusp in their density profile: while we search for the bestfitting King model to the entire density profile (including the central cusp), in the Harris catalog the concentration parameter is arbitrarily fixed to a constant value (c=2.5) and R c is determined as the radius where the surface brightness is equal to half the central value (see Trager et al. 1993).Hence, for the purposes of our analysis, the simulated clusters at 13 Gyr can be considered well representative of the Galactic GC population.

Normalised cumulative radial distributions and nCRD parameters
Once the value of r h is obtained from the King model fit to the projected density profile, to determine the new diagnostics of dynamical evolution defined in B22 and B23, the first step is to build the nCRDs of all the stars with V < V TO + 0.5, and located within a projected distance (R) equal to 0.5 × r h from the center, r h being the three-dimensional half-mass radius of the singlemass King model fitting the projected density profile (built by using stars with V < V T O + 1).
Figures 4 and 5 show the nCRDs of all the analyzed snapshots started with W 0 = 5 and W 0 = 7, respectively, with the right panels zooming into the innermost cluster region (R < 0.1r h ).The color code is based on the dynamical evolutionary stage identified for each extracted snapshot from the respective r 1% evolution (see caption).The dynamical evolution classification based on r 1% is well mirrored by different shapes of the nCRDs, with dynamically younger systems show-Bhat et al. ing systematically shallower curves, i.e., by construction, smaller percentages of stars included within the inner cluster regions.This is well in agreement with the expected effect of internal dynamical evolution, which tends to progressively increase the central density of star clusters, therefore steepening their nCRDs.The only exception (see the green line close to the red and blue ones in the right panel of Fig. 5) seems to be simulation N1 K0 RG6 TF0.025 W07 that is classified as pre-CC because, based on its r 1% evolution, shows no evidence of being close to CC, but has a nCRD "shifted within the group of highly evolved systems".A closer inspection of the figure, however, shows that this nCRD is the least steep within this group.This likely indicates that the system is more dynamically evolved than all the other ones classified as pre-CC, but still in an earlier phase than what we labelled "nCC".The three parameters presented in B22 and B23 have been specifically defined to quantify the observed morphological changes of the nCRD, which, in turn, are tracers of the dynamical age of the systems.Here we briefly recall their definition: • A 5 is the area subtended by each nCRD between the center and a projected distance equal to 5% r h (R/r h = 0.05); • P 5 , is the value of the nCRD (i.e., the fraction of stars) at the same distance from the center; • S 2.5 is the slope of the straight line tangent to the nCRD at R/r h = 0.025.More specifically, it is the slope of the tangent to the third-order polynomial function that best-fits the nCRD.
Adopting the above definitions, we determined the values of A 5 , P 5 , and S 2.5 for all the snapshots extracted at 13 Gyr from our Monte Carlo simulation survey.Out of 92 snapshots, two have less than ten stars inside their 5% r h , thus preventing any meaningful measure of the nCRD parameters.Hence, they have been discarded and only 90 snapshots have been used in the following analysis.

Reference boundaries and models from previous works
For the proper interpretation of the present results, we take advantage of our previous works.As discussed in B23, star clusters in different dynamical phases tend to occupy different regions of the A 5 −P 5 , A 5 −S 2.5 , and P 5 − S 2.5 diagrams, irrespective of the primordial binary fractions and dark remnant content.We therefore exploit those results to identify regions in these diagrams where dynamically-young and dynamically-old systems are expected to lie.
In Figure 6 we plot the values of the three parameters measured in B22 and B23 for simulations with 0%, 10%, 20% primordial binary fractions (BF0, BF10, and BF20, respectively), and for a case where a large amount of BHs was retained within the cluster for its entire evolution (DRr).Different colors correspond to different dynamical stages: green, cyan, blue and orange for early, intermediate, CC, and post-CC stages, respectively.Indigo is used for the DRr simulation that shows no evolution from the dynamical point of view, due to heating effects of the retained BHs.It is evident that the snapshots that are in early and intermediate dynamical stages, as well as the snapshots of the DRr simulation, all occupy the lower-left part of the plots.Conversely, the snapshots in CC and post-CC phases Bhat et al. .P5 and parameters plotted against A5 (left and right panels, respectively) for the reference simulations discussed in B22 and B23, started with 0%, 10% and 20% binary fractions (BF0, BF10, and BF20, respectively) and retaining a large population of dark remnants for the entire time evolution (DRr).The color code of the symbols is as follows: green, cyan, blue, and orange for early, intermediate, CC, and post-CC phases, respectively; indigo diamonds refer to the DRr simulation, which always remains in very early phases of dynamical evolution.The boundary region including only dynamically-young (pre-CC) systems is shaded in green.The polynomial fits to the points are plotted as black dashed lines.
occupy the upper-right part of the diagrams.There is also a small overlapping region where a few intermediate and highly evolved snapshots are found to lie together.Based on this evidence, we drew a boundary region (green shaded) safely including only dynamicallyyoung snapshots from all the simulations analyzed in previous works.Specifically, the boundaries for pre-CC systems are set at A 5 ≤ 0.0015, P 5 ≤ 0.085 and S 2.5 ≤ 1.9.Although these boundaries are somehow arbitrary, they serve as reference and guide for the interpretation of the results obtained from the present simulation survey.For the same purposes we also fitted a polynomial function to the sequences of A 5 versus P 5 , and A 5 versus S 2.5 (black dashed lines in Fig. 6).

Effects of different initial conditions
The dynamical evolution of star clusters is known to be affected by the combination of several different processes and cluster properties.Here, we take advantage of the simulation survey to search for possible dependences of the results on each initial condition parameter considered individually, namely, the total number of stars N , the BH kick velocity prescription (k), the galactocentric distance (R g ), the tidal filling parameter (TF), and the King dimensionless potential (W 0 ; see Table 1).
Figure 7 summarizes the results obtained for the A 5 parameter.These are perfectly analogous to those found for P 5 and S 2.5 , which are therefore not shown here.The dashed lines correspond to the value A 5 = 0.0015 that, based on our previous works, marks the separation between systems that are in early and in advanced dynamical stages.As a general consideration, the ma-jority of the snapshots are found in pre-CC dynamical phases (i.e., show A 5 ≤ 0.0015), while only a few (13 out of 90) are in the dynamically-old region (above the dashed line in the figure).This is consistent with the shapes of the nCRDs discussed in Section 3.3, most of which are in the "shallow group" (green lines).The top-left panel of Fig. 7 shows a mild indication that the most massive systems (N = 1M) do not reach the highest values of A 5 (corresponding to the highest levels of dynamical evolution), compared to simulations started with lower values of N .This indicates that, as expected, N plays an important role in the clusters' dynamical history; for a given initial spatial distribution, the halfmass relaxation time increases with the cluster mass, and thus the dynamical evolution proceeds on a longer timescale leading to the trend shown in this panel.A hint of more advanced dynamical stages for larger fractions of systems with small initial Galactocentric distances (R g = 2 and 4 kpc, compared to R g = 6 kpc) is observable in the top-right panel, while more significant trends are found between A 5 and the remaining parameters.Faster dynamical evolutions (i.e., larger values of A 5 and larger percentages of snapshots above the dashed lines) are observed for systems that starts with smaller BH retention (k = 0, compared to k = 1), higher initial compactness (i.e., smaller degree of tidal underfilling: TF=0.025, compared to TF=0.05 and 0.1), and higher concentration (W 0 = 7, compared to W 0 = 5).
In particular, while the snapshots with k = 0 cover the entire range of A 5 values, those with reduced BH kick velocity (k = 1) are essentially all clumped in the low-end of the distribution, corresponding to young dynamical .A5 parameter measured in the analyzed snapshots plotted as a function of the adopted initial conditions (see Table 1): from top-left, to bottom-right the total number of stars N , the BH kick velocity prescription (k), the galactocentric distance (Rg), the tidal filling parameter (TF), and the initial King dimensionless potential (W0).The dashed lines mark the value of A5 adopted as boundary to separate clusters in pre-CC and more advanced dynamical phases (A5 = 0.0015).
evolutionary states.This is due to the effects of retained BH populations, which act as energy source that powers core expansion (e.g., Heggie & Hut 2003), delaying the contraction of the inner region and extending the time needed to reach CC.The only two k = 1 snapshots that fall in dynamically-old region are those already discussed in Section 3 (N500 K1 RG2 TF0.025 W05 and N500 K1 RG2 TF0.025 W07).They show essentially the same value of A 5 (0.001995 and 0.001994, respectively) and their relatively fast dynamical evolution is due to the fact that all the BHs are ejected within the first few Gyr, thus leaving to the cluster inner regions enough time to substantially contract and reach CC within 13 Gyr.It is interesting to point out that values of A 5 corresponding to dynamically-young clusters (in a pre-CC phase) are not exclusively present in models retaining BHs (i.e., models with k = 1) but also in a large number of models with no BH retention (with k = 0).A strong impact on A 5 is also observed from the TF parameter: only simulations that started with TF= 0.025 (the most compact in the survey) reach the dynamically-old region in the plot, while for higher val-ues of TF, almost all the snapshots have A 5 ≤ 0.0015 and several clusters with TF=0.1 dissolve within 15 Gyr of evolution.

nCRD parameters as diagnostics of dynamical evolution
The adopted reference boundaries for pre-CC systems (green shaded regions in Figure 6) have been chosen based on snapshots extracted at different chronological times during the cluster evolution, from the simulations discussed in B22 and B23.Conversely, in the present work the A 5 , P 5 and S 2.5 parameters have been all measured in snapshots extracted at the same age (13 Gyr) from cluster simulations run from different initial conditions.Hence, it is not obvious a priori that the current results closely follow the previous ones.
To verify this possibility, Figure 8 shows the values of P 5 and S 2.5 as a function of A 5 (left and right panels, respectively) measured from the simulation survey discussed in this work, compared to the boundary region for pre-CC clusters defined from B22 and B23 (green shaded area) and the polynomial functions fitted to our previous results (dashed lines; see Sect.4.1).The present results  and 5, indicating different dynamical phases: green for dynamically-young snapshots (pre-CC), red and blue for systems in the verge to collapse (nCC) and that already reached CC, respectively.The snapshots with IMBH are encircled in black.The green shaded regions and the black dashed lines are, respectively, the reference boundaries for pre-CC systems and the polynomial fits to the simulations discussed in B22 and B23, as defined in Section 4.1 (the same as in Fig. 6).

Bhat et al.
(empty circles) are colored as in Figs. 4 and 5 based on the evolution of the 1% Lagrangian radius, with green color for the pre-CC snapshots, red color for the clusters close to CC, and blue for those that already reached CC at 13 Gyr.As apparent, they are all aligned along narrow sequences that closely follow the polynomial fits to the results of B22 and B23.In addition, the green shaded region is populated exclusively by snapshots classified as pre-CC, while all the dynamically-old systems lie outside this boundary.Also the simulations classified as nCC show values outside the reference boundary for pre-CC snapshots, indicating their higher level of dynamical evolution, which is also testified by the shape of their nCRDs, despite their 1% Lagrangian radius did not formally reach CC at 13 Gyr yet.Although a clear distinction among the different stages of high dynamical evolution (nCC, CC, post-CC) is not possible from the nCRD parameters, this figure clearly demonstrates the suitability of the new diagnostics and the adopted boundaries to substantially distinguish between dynamically-young and dynamically-old snapshots.
The green circle close to the top-right corner of the boundary region is simulation N1 K0 RG6 TF0.025 W07, the one with the nCRD "shifted within the group of highly evolved systems" (right panel of Fig. 5).Its position in this diagram suggests that this system is more evolved than all the other pre-CC clusters, fully confirming the conclusions drawn in Section 3.3.In turn, this demonstrates that the nCRD parameters have higher sensitivity to dynamical evolution than r 1% , which is unable to quantify any difference between this case and the other pre-CC ones.
Figure 9 shows A 5 plotted against the King concentration parameter (c) determined from the King model fitting to the density profile of each analyzed snapshot (colored circles).For reference, the small grey circles are the values of A 5 obtained by integrating the density profile of a sequence of King models with concentration parameter c varying between 1 and 2.5, in steps of 0.05.As apparent, the values obtained for the dynamically young snapshots (green circles) well follow the sequence expected for King models up to c≃ 1.8, with some scatter.Then, the values of A 5 measured in dynamically-old systems (red and blue circles) systematically and significantly exceed those expected for a King model with the same concentration c.This is in qualitative agreement with the well-known behavior of the density profile, which significantly deviates from the King model distribution (due to the presence of a central power-law cusp) in advanced stages of dynamical evolution.Very interestingly, also the pre-CC cluster that is close to the boundary between dynamically-young and dynamically-old systems (simulation N1 K0 RG6 TF0.025 W07) shows a significant deviation from the King model sequence at c > 1.8, further confirming that it is more dynamically evolved than its pre-CC peers.
The systems hosting a central IMBH tend to have larger King concentrations (between 1.63 and 1.85) than most of the other pre-CC clusters, but closely follow the King sequence.This indicates that their density profile shows no significant deviations with respect to the King model behavior.Indeed, previous works (e.g., Baumgardt et al. 2005;Miocchi 2007, but see also Vesperini & Trenti 2010) predict that a central IMBH should induce a cusp in the innermost portion of the density profile, with a slope that is much shallower than the one expected in CC systems.

SUMMARY AND DISCUSSION
In this paper we have measured the three nCRD parameters proposed by B22 and B23 (namely, A 5 , P 5 and S 2.5 ) in a sample of 90 GCs simulated through the MOCCA Monte Carlo code, starting from different initial conditions (see Table 1).The simulated clusters have been analyzed following the methods commonly implemented in observational studies.The crucial difference with respect to our previous studies, is that in B22 and B23 the parameters were measured for just a few models at different evolutionary times corresponding to different chronological ages (and therefore with snapshots including stars spanning different mass ranges), while here we consider a much broader range of initial conditions and compare all the systems at the same chronological age (13 Gyr).Thanks to the different intial conditions these systems have reached different dynamical evolutionary stages, as clearly testified by the variety of shapes observed for the time evolution of their 1% Lagrangian radii (see a few examples in Figs. 1  and 2), and the sample therefore is more appropriate to represent the Galactic GC population, where all systems are coeval and old, with ages of 10-13 Gyr, but are characterized by distinct dynamical ages (see, e.g., Ferraro et al. 2018).Indeed, the comparisons between the structural parameters (King concentration, core and half-light radii) and the core relaxation time of the synthetic and the observed clusters show that these simulations are in general well representative of the Galactic GC population.
We find that the three nCRD parameters measured in the present sample are fully consistent with those obtained in B22 and B23, drawing the same sequences in the P 5 and S 2.5 versus A 5 diagrams, and with the pre-CC systems showing systematically smaller val-ues than the more dynamically evolved clusters (Fig. 8).The dependence the nCRD parameters on the initial conditions of the simulations (see Fig. 7 for the case of A 5 ) shows that, as expected, clusters that form more massive and less compact (i.e., those having longer initial half-mass relaxation times) tend to experience a slower dynamical evolution.This is also the case for clusters retaining larger populations of BHs until the present time.Some compact clusters (specifically the models N500 K1 RG2 TF0.025 W05 and N500 K1 RG2 TF0.025 W07) retain a large fraction of BHs initially, but these are then dynamically ejected leaving few or no BHs at the present time; in those cases the early ejection of BHs leaves enough time for the cluster to progressively contract and approach CC under the effect of two-body relaxation (see also Weatherford et al. 2018;Kremer et al. 2020, B23).The present-day properties thus depend on the interplay between initial structural properties and their effects on the clusters' stellar content, and dynamical energy sources such as BHs and binary stars.The clusters classified as dynamically-old show large King concentration (c > 1.8) and nCRD parameters that are systematically larger than the values expected from the direct integration of the King models computed with the values of c providing the best fit to their density profiles (see Fig. 9).This is in agreement with the well-established presence of a power-law cusp (i.e., a significant deviation from the King model behavior) in the central portion of the density profile of highly evolved GCs.
A further demonstration of the effectiveness of the nCRD parameters as diagnostics of internal dynamical evolution is provided in Figure 10, which shows a clear anti-correlation between A 5 and the present-day value of the central relaxation time (t rc ).It is interesting to note that nCC clusters are better distinguishable from CC systems in this plot, compared to Figs. 8 and 9.In fact, for the same value of A 5 they show systematically larger values of t rc , consistent with their lower dynamical evolution.In addition, A 5 is able to distinguish nCC clusters from pre-CC snapshots (above and below the threshold, respectively) that, showing essentially the same values of t rc , would not be separable in based on their central relaxation time.Very interestingly, the snapshots hosting an IMBH in their center are all clustered in the upper-left end of the pre-CC sequence in Fig. 10, at 0.00045 < A 5 < 0.001 and log(t rc ) ∼ 7.5.The IMBHs have masses ranging between 0.07% and 0.45% of the total cluster mass (i.e., between ∼ 100 and 1000M ⊙ ), and they are bound in binary systems with a neutron star, a white dwarf or a non-degenerate star, but in all cases no other BH is present in the system.Instead, all the systems currently hosting large populations of stellar-mass BHs (more than ∼ 30) are found at much larger values of t rc and smaller values of A 5 (green circles filled in grey in Fig. 10).This suggests that the heating effect of a few hundred solar mass IMBH alone tends to be weaker than that of several stellar-mass BHs, and this result might help breaking the degeneracy between the presence of a single massive compact object and a concentration of dark remnants in the center of a GC, which often remain as possible indiscernible alternatives to explain an observed rising of the velocity dispersion in the innermost cluster region (see, e.g., van den Bosch et al. 2006 9 (also see the legend), with the addition of the green circles filled in grey that correspond to clusters hosting more than 30 stellar-mass BHs at an age of 13 Gyr.
All these findings prove that the three nCRD parameters A 5 , P 5 and S 2.5 are powerful diagnostics of the stage of internal dynamical evolution reached by star clusters.Compared to the computation of t rc , they offer the advantage of being fully empirical and not relying on simplifying assumptions and approximations.With respect to alternative methods requiring the determination of the stellar mass function or internal kinematics at different radial distances from the center (e.g., Baumgardt & Makino 2003;Tiongco et al. 2016;Bianchini et al. 2016;Bianchini et al. 2018;Webb & Vesperini 2017), the diagnostics presented here and those based on exotic populations as blue straggler stars (e.g., Ferraro et al. 2018Ferraro et al. , 2020) ) have the advantage to be more easily measurable from observational data.In fact, they just require high-resolution photometry down to ∼ 0.5 magnitudes below the main-sequence turnoff point and within a projected distance equal to 0.5 − 1 × r h from the center.This makes them suitable to recognize highly evolved systems also in cases where the central density cusp is not detectable from observations (e.g., Ferraro et al. 2023), either because of too scarce statistics in the construction of the density profile, or because the cluster is not CC yet or is in a rebound stage during its post-CC evolution.Because of their relatively simple determination, they also appear to be promising tools for the investigation of the dynamical age of stellar systems even beyond the Milky Way, out to distances where individual stars in the upper main-sequence can be resolved.A work specifically devoted to the observational determination of the nCRD parameters in a large sample of Galactic GCs is in preparation.

Figure 1 .Figure 2 .
Figure1.Time evolution of the 1% Lagrangian radius for the simulation runs N500 K1 RG2 TF0.025 W05, N750 K0 RG2 TF0.05 W07, N750 K0 RG4 TF0.1 W05 and N1 K0 RG2 TF0.1 W07 (from top left, to bottom right).The red vertical line at 13 Gyr marks the time snapshot of the simulation used for analysis.The classification of the 13 Gyr time snapshots in terms of dynamical evolutionary stage (CC, near to CC, and pre-CC) is marked in each panel.

Figure 3 .
Figure3.Half-light radius (R hl ) plotted against core radius (Rc, left panel), core radius plotted against the King concentration parameter (c, right panel) and core relaxation time (trc) plotted against concentration for the 92 simulation snapshots extracted at 13 Gyr (cyan squares for W0 = 5, and red circles for W0 = 7), compared to those measured in the Galactic GC population (fromHarris 1996), which are shown as filled circles (gray color for the CC and post-CC systems with c=2.5, black color for the others).

Figure 4 .Figure 5 .
Figure4.Left panel: Normalized cumulative radial distributions (nCRDs) of the stars brighter than VTO + 0.5 and located within 0.5×r h from the center for the snapshots extracted at 13 Gyr from the W0 = 5 runs.Right panel: zoom into the inner cluster region.The color code is as follows: green for dynamically young (pre-CC) snapshots, red for snapshots close to CC (nCC), and blue for dynamically-old snapshots.
Figure6.P5 and parameters plotted against A5 (left and right panels, respectively) for the reference simulations discussed in B22 and B23, started with 0%, 10% and 20% binary fractions (BF0, BF10, and BF20, respectively) and retaining a large population of dark remnants for the entire time evolution (DRr).The color code of the symbols is as follows: green, cyan, blue, and orange for early, intermediate, CC, and post-CC phases, respectively; indigo diamonds refer to the DRr simulation, which always remains in very early phases of dynamical evolution.The boundary region including only dynamically-young (pre-CC) systems is shaded in green.The polynomial fits to the points are plotted as black dashed lines.
Figure7.A5 parameter measured in the analyzed snapshots plotted as a function of the adopted initial conditions (see Table1): from top-left, to bottom-right the total number of stars N , the BH kick velocity prescription (k), the galactocentric distance (Rg), the tidal filling parameter (TF), and the initial King dimensionless potential (W0).The dashed lines mark the value of A5 adopted as boundary to separate clusters in pre-CC and more advanced dynamical phases (A5 = 0.0015).

Figure 8 .
Figure8.P5 and S2.5 parameters plotted against (left and right panels, respectively), as measured in the 13 Gyr snapshots of the simulation survey discussed in this work.The color code is the same as in Figs.4 and 5, indicating different dynamical phases: green for dynamically-young snapshots (pre-CC), red and blue for systems in the verge to collapse (nCC) and that already reached CC, respectively.The snapshots with IMBH are encircled in black.The green shaded regions and the black dashed lines are, respectively, the reference boundaries for pre-CC systems and the polynomial fits to the simulations discussed in B22 and B23, as defined in Section 4.1 (the same as in Fig.6).

Figure 9 .
Figure9.A5 measured in the survey simulation snapshots (large circles) plotted against their concentration parameter as obtained from the best-fit King model to their density profile.The color code is the same as in Fig.8, and the green dashed line marks the separation boundary between pre-CC and more evolved clusters.The small grey circles are the values of A5 obtained by directly integrating the density profile of King models with concentration c varying in the range 1 and 2.5, in steps of 0.05.The bottom panel shows the residuals between the measured nCRD parameter and the value expected for the King model with same concentration.
for the case of M15 and Vitral et al. 2023 for M4).

Figure 10 .
Figure10.Logarithm of A5 plotted against the logarithm of trc.The symbols are as in Fig.9(also see the legend), with the addition of the green circles filled in grey that correspond to clusters hosting more than 30 stellar-mass BHs at an age of 13 Gyr.

Table 1 .
Initial conditions of the simulations