EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 698 (June 2025) A&A, 698 (2025) L11 Full HTML
Open Access
Issue
A&A
Volume 698, June 2025
Article Number L11
Number of page(s) 7
Section Letters to the Editor
DOI https://doi.org/10.1051/0004-6361/202554233
Published online 03 June 2025

© The Authors 2025

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1. Introduction

Globular clusters (GCs) were historically thought to host simple stellar populations with uniform chemical compositions. However, photometric and spectroscopic studies over the past few decades have revealed the presence of multiple stellar populations (MPOP), primarily characterized by internal variations in light-element abundances (see Bastian & Lardo 2018; Gratton et al. 2019, and references therein). These populations include a pristine first population (1P), whose chemical composition resembles that of field stars with similar metallicities, and one or more enriched populations (2P), which exhibit signatures of high-temperature hydrogen burning (Milone & Marino 2022). The fraction of enriched 2P stars in Milky Way (MW) GCs strongly correlates with the present-day cluster mass, ranging from approximately 40% in low-mass clusters to nearly 90% in the most massive ones (see Milone & Marino 2022, and references therein).

The formation of MPOP in GCs remains a topic of active debate, with several theoretical scenarios proposed to explain its origin. However, none of these scenarios fully accounts for all the observational features reported in the literature. Giersz et al. (2025) provide a concise overview of the most widely accepted formation mechanisms, particularly from the perspective of numerical simulations. Their study also suggests that the post-formation migration of GCs may play a crucial role in shaping the properties of MPOP. Most formation scenarios, including those that attribute MPOP to pollution from asymptotic giant branch (AGB) stars, predict that 2P should be more centrally concentrated than 1P. Observations generally support this prediction, but exceptions exist, such as the Galactic GCs NGC 3201 and NGC 6101, where 1P red giant branch (RGB) stars are found to be more centrally concentrated than 2P RGB stars (Leitinger et al. 2023). However, in the case of NGC 3201, Mehta et al. (2024), using RGB stars along with other evolved stars from the horizontal branch (HB), found no evidence of such overconcentration, as 1P and 2P appear to be well mixed.

In this short paper, we present results from a GC model simulated with the MOCCA code, in which the 2P population initially has a much higher central concentration than the 1P. However, by 13 Gyr, the 1P RGB stars become more spatially overconcentrated than the 2P RGB stars. The overconcentration of 1P RGB stars relative to 2P RGB stars may provide valuable insights into their formation scenarios and long-term evolution. It could also serve as an important constraint on the initial conditions of GC formation. Therefore, understanding the mechanisms driving this phenomenon and identifying the clusters where it may occur are crucial. In this study, we focus on the transient nature of 1P RGB overconcentration, offering a general explanation for its formation and the types of GCs in which it may be observed. A more detailed investigation will be presented in a forthcoming paper.

The paper is structured as follows. Section 2 provides a brief overview of the MOCCA code. Section 3 presents the simulation results for models exhibiting 1P overconcentration, and Section 4 interprets these findings and discusses their implications.

2. Method

This study is based on simulations carried out with the MOCCA Monte Carlo code, which models the stellar and dynamical evolution of realistic GCs, including multiple populations and delayed 2P formation. For details of the code and physical assumptions, see Appendix A.

The model presented in this paper has the following initial conditions: the number of objects, N1P = 800 000, N2P = 300 000, and a Kroupa (2001) initial mass function (IMF), with 1P stars sampled between 0.08 M⊙ and 150 M⊙, and 2P stars between 0.08 M⊙ and 20 M⊙. The metallicities are Z1P = 0.001 and Z2P = 0.001, with binary fractions fb1P = 0.95 and fb2P = 0.951. The Galactocentric distance is Rg = 1.0 kpc, and the ratio of half-mass radii is Rh2P/Rh1P = 0.1. 1P and 2P are separately in virial equilibrium, with 1P being tidally filling (with Rh1P = 9.03 pc and Rt/Rh1P = 3.66). The initial King (1966) model concentration parameter is W0 = 3 for 1P and 7 for 2P. Over 1 Gyr following 5.5 Gyr of evolution, the cluster migrates along its circular orbit from 1 kpc to 5 kpc, mimicking the capture of a GC during a dwarf galaxy merger with the MW.

3. Results

Inspired by the work of Leitinger et al. (2023), who conducted a homogeneous analysis of MPOP in 28 Galactic GCs using a combination of HST photometry and wide-field, ground-based photometry of RGB stars, we examine a puzzling trend that they identified. Their study found that in a few Galactic GCs, such as NGC 3201 and NGC 6101, 1P RGB stars are more concentrated within a few times the half-light radius than 2P RGB stars. This unexpected result is difficult to reconcile with current MPOP formation scenarios. Recently, Mehta et al. (2024) cast doubt on these findings, showing that in NGC 3201, when HB and AGB stars are included alongside RGB stars, both populations exhibit similar spatial distributions. However, Cadelano et al. (2024) confirmed the results of Leitinger et al. (2023) for NGC 3201, though they found that 2P RGB stars follow a bimodal distribution and are more concentrated within the half-light radius.

With these results in mind, we took a closer look at the MPOP simulations performed with the MOCCA code, as is discussed in Giersz et al. (2025). We identified several models in which 1P stars are more concentrated than 2P stars. Studies by Geisler et al. (1995), Massari et al. (2019), Chen & Gnedin (2024) suggest that NGC 6101 and NGC 3201 likely formed ex situ, a scenario supported by the fact that both clusters follow retrograde orbits. Based on this, we selected a GC model with a capture time by the MW of 5.5 Gyr and an orbital migration time of 1 Gyr. In addition, we tried to select a model with observational properties similar to those of NGC 3201, specifically: a cluster mass of 1.45 × 105 (1.93 × 105) M⊙, a core radius of Rc = 1.22 (1.66) pc, an observational core radius of Rc, obs = 3.13 (1.85) pc, a half-light radius of Rhl = 6.22 (5.17) pc, a half-mass radius of Rh = 8.71 (8.68) pc, a tidal radius of Rt = 81.72 (82.34) pc, and the ratio for RGB stars between the 1P and 2P populations of 0.59 (0.46). The values for the MOCCA model were taken at 13 Gyr, while the values in parentheses represent the observed properties of NGC 3201 from Baumgardt & Hilker (2018) and Milone et al. (2017). The observational properties of the NGC 3201 are close to those of the model (except for Rc, obs – see Appendix C for a discussion), but this does not mean that the model represents the actual evolution of the GC. It only serves as a possible example of GC evolution, and we should not expect all cluster properties to be the same, in particular, the 1P concentration range. The model is used to illustrate the possible causes of the observed 1P and 2P properties in NGC 3201.

Figure 1 shows the cumulative number distributions of RGB stars for 1P and 2P as a function of distance scaled by Rhl at several selected snapshot times for a chosen MOCCA model. It clearly shows that 1P is more centrally concentrated than 2P only in the snapshot at 13 Gyr. In other snapshots, both populations are well mixed, or 2P is more concentrated than 1P. The number distributions of both populations vary significantly over time, indicating that the overconcentration of 1P is a transient feature. The magnitude of 1P overconcentration at 13 Gyr is not as pronounced as that observed for NGC 3201 by Leitinger et al. (2023) and Cadelano et al. (2024). This is not unexpected, as the presented model is only a rough approximation of an actual GC. As is discussed in Section 4, the extent of 1P overconcentration depends on the evolutionary details of RGB progenitors, the migration time, and the initial conditions of the cluster model.

thumbnail Fig. 1.

Cumulative number distributions of RGB stars for 1P (black line) and 2P (red line) as a function of distance scaled by Rhl for selected projected (2D) snapshots from 12 Gyr to 13.5 Gyr. 1P is more centrally concentrated than 2P only in the snapshot at 13 Gyr. The number distributions of both populations change significantly over time, and the overconcentration of 1P is a transient feature. The Rhl values, snapshot times, and the number of RGB stars are provided in the insets of each panel.

Since the overconcentration of 1P relative to 2P appears to be a transient phenomenon, we examined whether this effect is also present when considering stellar types other than RGB stars. We selected the snapshot at 13 Gyr, and as in Figure 1, all distances are scaled by Rhl. Figure 2 presents cumulative distributions for different stellar populations. We analyzed main sequence (MS) stars with masses greater than 0.4 M⊙ (a mass range typically observable in most GCs) as well as all stars, including compact remnants, belonging to 1P and 2P. The top left panel shows the cumulative mass distribution of all 1P and 2P stars, while the top right panel presents the cumulative number distributions of all stars for both populations. The bottom left panel displays the cumulative number distribution of MS stars with masses greater than 0.4 M⊙ for 1P and 2P. Finally, the bottom right panel presents the cumulative number distributions of RGB stars for 1P and 2P up to a distance equal to Rhl. The overconcentration of 1P relative to 2P is also evident when analyzing the mass distribution of all stars in the cluster. While the overconcentration is subtle, it remains noticeable. For both populations, the half-mass radii were already mixed at 5 Gyr. However, when considering the number distribution of stars instead of the mass distribution, the two populations appear to be perfectly mixed. From the analysis of the model, it appears that stellar-mass black holes (BHs) are responsible for this effect. These BHs originate from 1P stars, as 2P stars have a maximum initial mass of only 20 M⊙, which significantly limits 2P BH formation due to the IMF. Over time, these BHs form a subsystem (BHS) at the cluster center. The influence of the BHS on the global evolution of the cluster is discussed in Section 4.1. When considering only MS stars that are observable in distant GCs (e.g., those with masses greater than 0.4 M⊙), the overconcentration of 1P stars in the number distribution is not apparent, as both populations appear well mixed. This suggests that the type and number of selected stars play a crucial role in detecting the overconcentration of 1P relative to 2P. Leitinger et al. (2023) and Cadelano et al. (2024) reported that in the central region of NGC 3201 (inside Rhl), 2P stars are more concentrated than 1P stars or are nearly mixed. A similar trend is observed in the bottom right panel of Figure 2. In the next section (Section 4), we further analyze the properties of models that explain why 1P is more concentrated than 2P, why the overconcentration of 1P is observed only in RGB stars, and why this feature appears to be transient.

thumbnail Fig. 2.

Top left panel: Cumulative mass distributions of all stars for 1P (black line) and 2P (red line) as a function of distance scaled by Rhl in the projected (2D) snapshot at 13 Gyr. Top right panel: Cumulative number distributions of all stars for 1P (black line) and 2P (red line) as a function of distance scaled by Rhl in the snapshot at 13 Gyr. Bottom left panel: Cumulative number distributions of MS stars with masses greater than 0.4 M⊙ for 1P (black line) and 2P (red line) as a function of distance scaled by Rhl in the snapshot at 13 Gyr. Bottom right panel: Cumulative number distributions of RGB stars for 1P (black line) and 2P (red line) as a function of distance scaled by Rhl in the snapshot at 13 Gyr, shown up to a distance of Rhl.

4. Discussion

4.1. Black hole subsystem and cluster migration

By analyzing MOCCA simulations with the MPOP implementation, as is described in Giersz et al. (2025), and focusing on models in which 1P is more centrally concentrated than 2P, we draw the following conclusions. The overconcentration of 1P is observed in GC models that have lost 80–90% of their initial mass, leading to present-day (11–13 Gyr) masses of approximately 105 M⊙. In all models that show 1P overconcentration, the feature is transient. Additionally, these clusters retain a BHS comprising up to a few hundred stellar-mass BHs at the present day. In clusters for which the present-day mass is only 10–20% of the initial mass and the BHS accounts for approximately 5% of the total cluster mass, the mass distribution of 1P becomes more centrally concentrated than that of 2P. All models with 1P overconcentration are characterized by a very wide range of tidal radii from 10–20 to about 80–90 pc. Only models with migration at a late time show properties similar to the ones observed in MW GCs.

The GC model described in this paper initially retains nearly 1000 stellar-mass BHs at around 50 Myr, originating from the evolution of massive 1P progenitors. At this time, the 10% and 50% Lagrangian radii for these BHs are approximately 4.8 pc and 10.8 pc, respectively, and begin to segregate toward the center, finally forming a BHS (see Appendix B for a more detailed description of the BHS evolution). The BHS, consisting of both single and binary BHs, strongly influences the cluster’s dynamical evolution through interactions with other stars (Mackey et al. 2008; Breen & Heggie 2013; Sippel & Hurley 2013; Morscher et al. 2015; Arca Sedda et al. 2018; Kremer et al. 2019; Weatherford et al. 2020). Since the 2P population is initially more centrally concentrated in this model, 2P stars are more likely to undergo strong encounters with 1P BHs and other 2P stars. These interactions can scatter 2P stars onto wider orbits or eject them from the cluster entirely (see Section 4.2).

For both NGC 3201 and NGC 6101, where 1P RGB stars appear more centrally concentrated than 2P RGB stars, several studies suggest that these GCs likely harbor a significant population of stellar-mass BHs (Peuten et al. 2016; Askar et al. 2018; Kremer et al. 2018; Vitral et al. 2022). Furthermore, three stellar-mass BHs have been identified in NGC 3201 through radial velocity variations using MUSE spectroscopic observations (Giesers et al. 2018, 2019). Therefore, it is possible that the possible presence of a BHS has played a role in shaping the spatial distribution of 1P and 2P RGB stars in these GCs.

The GC NGC 3201 is likely an ex situ cluster that was accreted by the MW along with a dwarf galaxy in the Gaia-Enceladus–Sequoia event (Massari et al. 2019). It is characterized by large Rh and Rt radii. During the early stages of its evolution, it likely experienced a strong tidal field within its parent galaxy, and it can be assumed that it was initially tidally filling, with small Rh and Rt radii (Renaud et al. 2017). After being captured by the MW, the tidal field weakened considerably, causing the cluster to become tidally underfilled, allowing it to expand freely within its Rt and significantly increase its half-mass relaxation time. For these reasons, we selected a cluster model that initially forms near the galactic center with an initial mass of approximately 106 M⊙. Over a few gigayears, the cluster migrates to a larger Galactocentric distance (see Section 3 for details). This model not only reproduces similar observational properties to NGC 3201 but also yields a comparable ratio of 2P to 1P stars. We would like to strongly emphasize that GC migration time is a free parameter in simulations and plays an important role in creating conditions for achieving 1P overconcentration. At the same time, the magnitude of 1P overconcentration depends on many initial parameters such as the number ratio of 2P to 1P, the mass of the GC, the overconcentration of 2P relative to 1P, the galactocentric position of the GC, and the strength of the galactic tidal field. The capture time of NGC 3201 by the MW is about 10 Gyr (Massari et al. 2019). This is about twice as large as is assumed in the model. Since the strength of the 1P overconcentration depends on many initial parameters, we are confident that the appropriate selection of these parameters will result in a model of the NGC 3201 cluster that corresponds much better to observations. This will be the subject of the next paper.

Due to the cluster’s relatively low present-day mass, the number of RGB stars is relatively small, on the order of several hundred. As we show below, such a small number of RGB stars plays a crucial role in the overconcentration of 1P relative to 2P RGB stars.

4.2. Preferential ejection of enriched RGB progenitors by strong dynamical encounters

To explain the spatial overconcentration of 1P RGB stars at 13 Gyr, we first identified stars that evolve into RGBs at this epoch. These stars are near the turn-off mass, which is 0.8154 M⊙ for Z = 0.001 and had a zero-age MS mass of approximately 0.83 M⊙. Assuming that all stars in the initial model (or at ∼100 Myr, when 2P is fully added) with masses between 0.81 and 0.83 M⊙ will have evolved into RGB stars by 13 Gyr, we defined this mass range as the present-day RGB progenitors and examined those that escaped before 13 Gyr.

For 1P stars, we counted all escaping RGB progenitors within the 0.81–0.83 M⊙ mass range, regardless of the cause of escape, and found that 5048 1P RGB progenitors had escaped the cluster. Similarly, for 2P stars, 1167 RGB progenitors had escaped. This significant difference is expected, as 1P was initially tidally filling. However, when considering only escapers due to strong dynamical interactions (binary-single or binary-binary encounters), the numbers shift dramatically: 1P RGB progenitor escapers drop to 159, while 2P increases to 414. Thus, 2P RGB progenitors are 2.6 times more likely to be dynamically ejected than 1P progenitors, which is consistent with their initial central concentration.

We find that 2P RGB progenitors are more frequently ejected from the inner regions of the cluster via strong dynamical encounters. As is shown in Fig. 3, most 2P escapers are removed between 2 and 6 Gyr, primarily from within a few times the half-light radius (Rhl). After 6 Gyr, following the cluster’s migration to a larger Galactocentric distance, the ejection rate declines but remains higher for 2P stars. The lower panel shows that 2P escapers are mostly launched from the inner regions, while the top histogram highlights their higher escape rate compared to 1P. Markers in the figure also reveal that a significant fraction of these ejections – especially for 2P stars – involve interactions with BHs (23.2% for 2P vs. 12.6% for 1P). Together, these results confirm that strong dynamical encounters, often involving BHs, preferentially remove 2P RGB progenitors, contributing to the overconcentration of 1P RGB stars at 13 Gyr.

thumbnail Fig. 3.

The lower panel: Escape time versus radial position for dynamically ejected present-day RGB progenitor stars with initial masses between 0.81 and 0.83 M⊙. Each point represents a present-day RGB progenitor that escaped due to a strong three-body (binary-single) or four-body (binary-binary) encounter. The x axis shows the escape time, and the y axis indicates the radial position within the cluster at which the escape-causing interaction occurred. Superimposed lines trace the evolution of the core radius (Rc, grey), half-mass radius (Rh, dark blue), half-light radius (Rhl, light blue), and tidal radius (Rt, brown). Violet points represent 1P escapers, and green points represent 2P escapers; filled symbols denote single stars, and open symbols denote binaries. Triangles indicate stars that escaped directly due to an encounter involving at least one BH, while squares indicate stars that had at least one such BH encounter during their dynamical history. Of the 159 1P escapers, 20 (12.6%) experienced BH interactions; among the 414 2P escapers, 96 (23.2%) were involved in BH interactions. The vertical dashed red line marks the 13 Gyr snapshot, where the 1P overconcentration is most evident. The top panel: Histograms of escape times for both populations, binned in 200 Myr intervals, highlighting the higher escape rate of 2P RGB progenitors from strong encounters.

To quantify this effect, we defined the ratio:

R = N RGB , in _ cluster N RGB , progenitor _ dyn _ escapers , $$ \begin{aligned} R = \frac{N_{\mathrm{RGB,in}\_{\rm cluster}}}{N_{\mathrm{RGB,progenitor}\_{\rm dyn}\_{\rm escapers}}}, \end{aligned} $$

where N RGB , in _ cluster $ N_{\mathrm{RGB,in}\_{\rm cluster}} $ is the number of RGB stars present in the cluster at a given snapshot, and N RGB , progenitor _ dyn _ escapers $ N_{\mathrm{RGB,progenitor}\_{\rm dyn}\_{\rm escapers}} $ is the cumulative number of RGB progenitors lost to strong encounters by that time. For 2P at 13 Gyr, we find R = 1.04 (430/414), indicating that the number of remaining RGB stars is nearly equal to those lost through strong encounters. At earlier times, such as 12.5 Gyr, this ratio is higher (R = 1.10), as fewer progenitors have escaped relative to those still in the cluster. By 13.5 Gyr, the ratio rises to R = 1.11, showing that the overconcentration of 1P relative to 2P RGB stars is most pronounced at 13 Gyr, when R is closest to unity. This near-parity magnifies the impact of missing 2P RGB progenitors on the observed spatial distributions. At 13 Gyr, the number of 2P RGB stars in the cluster (430) is nearly matched by the number of progenitors ejected via strong encounters (414). By contrast, at 12.5 Gyr, the higher ratio (R = 1.1) means more 2P RGB stars remain in the cluster, lessening this effect. By 13.5 Gyr, R increases again to 1.11, diminishing the overconcentration. This increase occurs because lower-mass RGB progenitors, more numerous due to the IMF and still present in the cluster, begin evolving onto the RGB at this stage, thereby increasing the numerator in the ratio.

This analysis highlights the distinct dynamical evolution of 2P RGB progenitors, which are preferentially ejected through strong interactions, particularly from the cluster’s inner regions. The interplay of small-number statistics and dynamical encounters thus plays a crucial role in shaping the observed RGB spatial distribution.

4.3. Conclusions

Our study suggests that the spatial distribution and, potentially, the kinematic properties of MPOP may depend on the type of stars observed. This effect is particularly relevant for GCs with present-day masses of a few 105 M⊙, which have retained only 10–20% of their initial mass. Despite the significant mass loss, such clusters may appear dynamically young due to BHS heating and migration to larger Galactocentric distances. When observations are limited to RGB stars, small-number statistics and dynamical interactions can distort their spatial distribution, leading to biased conclusions about the overall distribution of MPOP. The overconcentration of 1P RGB stars relative to 2P RGB stars appears to be a transient effect.

Confirming the findings of Leitinger et al. (2023) and Cadelano et al. (2024), that in some MW GCs, 1P RGB stars are more centrally concentrated than 2P RGB stars, requires observations of MS stars. Our MOCCA simulations show no 1P overconcentration relative to 2P among MS stars. Verifying this observationally would be a crucial step toward understanding MPOP properties and supporting the AGB scenario, and possibly other scenarios, of MPOP formation. Additionally, our results align with the kinematic analysis by Cadelano et al. (2024), which suggests that 2P stars initially formed more centrally concentrated than 1P stars, influencing their present-day spatial and kinematic distributions.

We have shown that the observed overconcentration of 1P relative to 2P can be explained within the AGB scenario, provided certain internal and external conditions are met. This effect appears to be transient, driven by fluctuations in the number and spatial distribution of RGB stars. Further analysis is needed to confirm these results, particularly regarding binary fraction distributions and velocity anisotropy in 1P and 2P. In an upcoming study, we shall conduct a more detailed comparison of MOCCA model results with observations from Leitinger et al. (2025) and Cadelano et al. (2024) to further investigate these trends.

5. Data availability

The data for the MOCCA model described in this paper, including the 13 Gyr snapshot, is available at https://doi.org/10.5281/zenodo.14908094. Additional data can be provided upon request to the corresponding authors.

1

Initial binary parameters were distributed according to Belloni et al. (2017).

Acknowledgments

We are grateful to the anonymous reviewer for their insightful comments, which helped improve the manuscript. MG, AH, GW, LH were supported by the Polish National Science Center (NCN) through the grant 2021/41/B/ST9/01191. AA acknowledges support for this paper from project No. 2021/43/P/ST9/03167 co-funded by the Polish National Science Center (NCN) and the European Union Framework Programme for Research and Innovation Horizon 2020 under the Marie Skłodowska-Curie grant agreement No. 945339. AH acknowledges support by the IDUB grant 140/04/POB4/0006 at Adam Mickiewicz University in Poznan, Poland. For the purpose of Open Access, the authors have applied for a CC-BY public copyright license to any Author Accepted Manuscript (AAM) version arising from this submission.

References

  1. Arca Sedda, M., Askar, A., & Giersz, M. 2018, MNRAS, 479, 4652 [NASA ADS] [CrossRef] [Google Scholar]
  2. Askar, A., Arca Sedda, M., & Giersz, M. 2018, MNRAS, 478, 1844 [NASA ADS] [CrossRef] [Google Scholar]
  3. Baker, J. G., Boggs, W. D., Centrella, J., et al. 2008, ApJ, 682, L29 [NASA ADS] [CrossRef] [Google Scholar]
  4. Banerjee, S., Belczynski, K., Fryer, C. L., et al. 2020, A&A, 639, A41 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  5. Bastian, N., & Lardo, C. 2018, ARA&A, 56, 83 [Google Scholar]
  6. Baumgardt, H., & Hilker, M. 2018, MNRAS, 478, 1520 [Google Scholar]
  7. Baumgardt, H., Hilker, M., Sollima, A., & Bellini, A. 2019, MNRAS, 482, 5138 [Google Scholar]
  8. Belczynski, K., Kalogera, V., & Bulik, T. 2002, ApJ, 572, 407 [NASA ADS] [CrossRef] [Google Scholar]
  9. Belczynski, K., Heger, A., Gladysz, W., et al. 2016, A&A, 594, A97 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. Belloni, D., Askar, A., Giersz, M., Kroupa, P., & Rocha-Pinto, H. J. 2017, MNRAS, 471, 2812 [NASA ADS] [CrossRef] [Google Scholar]
  11. Belloni, D., Kroupa, P., Rocha-Pinto, H. J., & Giersz, M. 2018, MNRAS, 474, 3740 [CrossRef] [Google Scholar]
  12. Breen, P. G., & Heggie, D. C. 2013, MNRAS, 432, 2779 [NASA ADS] [CrossRef] [Google Scholar]
  13. Cadelano, M., Dalessandro, E., & Vesperini, E. 2024, A&A, 685, A158 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  14. Chen, Y., & Gnedin, O. Y. 2024, Open J. Astrophys., 7, 23 [NASA ADS] [Google Scholar]
  15. Dickson, N., Smith, P. J., Hénault-Brunet, V., Gieles, M., & Baumgardt, H. 2024, MNRAS, 529, 331 [NASA ADS] [CrossRef] [Google Scholar]
  16. Fregeau, J. M., & Rasio, F. A. 2007, ApJ, 658, 1047 [Google Scholar]
  17. Fregeau, J. M., Cheung, P., Portegies Zwart, S. F., & Rasio, F. A. 2004, MNRAS, 352, 1 [CrossRef] [Google Scholar]
  18. Fryer, C. L., Belczynski, K., Wiktorowicz, G., et al. 2012, ApJ, 749, 91 [Google Scholar]
  19. Fuller, J., & Ma, L. 2019, ApJ, 881, L1 [Google Scholar]
  20. Geisler, D., Piatti, A. E., Claria, J. J., & Minniti, D. 1995, AJ, 109, 605 [NASA ADS] [CrossRef] [Google Scholar]
  21. Giersz, M. 1998, MNRAS, 298, 1239 [NASA ADS] [CrossRef] [Google Scholar]
  22. Giersz, M., Heggie, D. C., Hurley, J. R., & Hypki, A. 2013, MNRAS, 431, 2184 [NASA ADS] [CrossRef] [Google Scholar]
  23. Giersz, M., Askar, A., Hypki, A., et al. 2025, A&A, in press, https://doi.org/10.1051/0004-6361/202452945 [Google Scholar]
  24. Giesers, B., Dreizler, S., Husser, T.-O., et al. 2018, MNRAS, 475, L15 [NASA ADS] [CrossRef] [Google Scholar]
  25. Giesers, B., Kamann, S., Dreizler, S., et al. 2019, A&A, 632, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Gratton, R., Bragaglia, A., Carretta, E., et al. 2019, A&ARv, 27, 8 [Google Scholar]
  27. Harris, W. E. 1996, AJ, 112, 1487 [Google Scholar]
  28. Hobbs, G., Lorimer, D. R., Lyne, A. G., & Kramer, M. 2005, MNRAS, 360, 974 [Google Scholar]
  29. Hurley, J. R., Pols, O. R., & Tout, C. A. 2000, MNRAS, 315, 543 [Google Scholar]
  30. Hurley, J. R., Tout, C. A., & Pols, O. R. 2002, MNRAS, 329, 897 [Google Scholar]
  31. Hypki, A., & Giersz, M. 2013, MNRAS, 429, 1221 [NASA ADS] [CrossRef] [Google Scholar]
  32. Hypki, A., Giersz, M., Hong, J., et al. 2022, MNRAS, 517, 4768 [NASA ADS] [CrossRef] [Google Scholar]
  33. Hypki, A., Vesperini, E., Giersz, M., et al. 2025, A&A, 693, A41 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  34. Kamlah, A. W. H., Leveque, A., Spurzem, R., et al. 2022, MNRAS, 511, 4060 [NASA ADS] [CrossRef] [Google Scholar]
  35. King, I. R. 1966, AJ, 71, 64 [Google Scholar]
  36. Kremer, K., Ye, C. S., Chatterjee, S., Rodriguez, C. L., & Rasio, F. A. 2018, ApJ, 855, L15 [NASA ADS] [CrossRef] [Google Scholar]
  37. Kremer, K., Chatterjee, S., Ye, C. S., Rodriguez, C. L., & Rasio, F. A. 2019, ApJ, 871, 38 [CrossRef] [Google Scholar]
  38. Kroupa, P. 2001, MNRAS, 322, 231 [NASA ADS] [CrossRef] [Google Scholar]
  39. Leitinger, E., Baumgardt, H., Cabrera-Ziri, I., Hilker, M., & Pancino, E. 2023, MNRAS, 520, 1456 [NASA ADS] [CrossRef] [Google Scholar]
  40. Leitinger, E. I., Baumgardt, H., Cabrera-Ziri, I., et al. 2025, A&A, 694, A184 [Google Scholar]
  41. Mackey, A. D., Wilkinson, M. I., Davies, M. B., & Gilmore, G. F. 2008, MNRAS, 386, 65 [NASA ADS] [CrossRef] [Google Scholar]
  42. Massari, D., Koppelman, H. H., & Helmi, A. 2019, A&A, 630, L4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  43. Mehta, V. J., Milone, A. P., Casagrande, L., et al. 2024, MNRAS, submitted [arXiv:2406.02755] [Google Scholar]
  44. Milone, A. P., & Marino, A. F. 2022, Universe, 8, 359 [NASA ADS] [CrossRef] [Google Scholar]
  45. Milone, A. P., Piotto, G., Renzini, A., et al. 2017, MNRAS, 464, 3636 [Google Scholar]
  46. Morawski, J., Giersz, M., Askar, A., & Belczynski, K. 2018, MNRAS, 481, 2168 [Google Scholar]
  47. Morscher, M., Pattabiraman, B., Rodriguez, C., Rasio, F. A., & Umbreit, S. 2015, ApJ, 800, 9 [NASA ADS] [CrossRef] [Google Scholar]
  48. Peuten, M., Zocchi, A., Gieles, M., Gualandris, A., & Hénault-Brunet, V. 2016, MNRAS, 462, 2333 [CrossRef] [Google Scholar]
  49. Renaud, F., Agertz, O., & Gieles, M. 2017, MNRAS, 465, 3622 [NASA ADS] [CrossRef] [Google Scholar]
  50. Sippel, A. C., & Hurley, J. R. 2013, MNRAS, 430, L30 [NASA ADS] [Google Scholar]
  51. Tanikawa, A., Yoshida, T., Kinugawa, T., Takahashi, K., & Umeda, H. 2020, MNRAS, 495, 4170 [NASA ADS] [CrossRef] [Google Scholar]
  52. Vitral, E., Kremer, K., Libralato, M., Mamon, G. A., & Bellini, A. 2022, MNRAS, 514, 806 [NASA ADS] [CrossRef] [Google Scholar]
  53. Weatherford, N. C., Chatterjee, S., Kremer, K., & Rasio, F. A. 2020, ApJ, 898, 162 [NASA ADS] [CrossRef] [Google Scholar]
  54. Zocchi, A., Gieles, M., & Hénault-Brunet, V. 2019, MNRAS, 482, 4713 [Google Scholar]

Appendix A: MOCCA code

This work is based on the numerical simulations performed with the MOCCA Monte Carlo code (Giersz 1998; Hypki & Giersz 2013; Giersz et al. 2013; Hypki et al. 2022, 2025; Giersz et al. 2025). MOCCA is an advanced code that performs full stellar and dynamical evolution of real-size star clusters up to the Hubble time. The implementation of stellar and binary evolution within the MOCCA code is based on the rapid population synthesis code BSE code (Hurley et al. 2000, 2002) that has been strongly updated by Belloni et al. (2017, 2018), Banerjee et al. (2020), and Kamlah et al. (2022). Strong dynamical interactions in MOCCA are performed with FEWBODY code (Fregeau et al. 2004; Fregeau & Rasio 2007). MOCCA can follow the full dynamical and stellar evolution of MPOP, allowing for delayed 2P formation relative to 1P. In this scenario, 2P forms through the re-accretion of gas surrounding the newly formed cluster, which mixes with enriched material from the ejected envelopes of AGB stars. Additionally, MOCCA incorporates the effects of cluster migration, driven by the rapidly evolving gravitational potential of the early Galaxy and dynamical friction.

The stellar and binary evolution parameters used in the presented simulations are referred to as Model C in Kamlah et al. (2022). In short, the metallicity of both populations in all the simulated models was set to Z = 0.001. The updated treatment for the evolution of massive stars was used according to Tanikawa et al. (2020) together with improved treatment for mass loss due to stellar winds and the inclusion of pair and pulsational pair-instability supernova (Belczynski et al. 2016). The masses of BHs and NS were determined according to the rapid supernovae prescriptions from Fryer et al. (2012). NS natal kicks were sampled from a Maxwelllian distribution with σ = 265 km/s (Hobbs et al. 2005). However, for BHs, these natal kicks were reduced according to the mass fallback prescription (Belczynski et al. 2002; Fryer et al. 2012). The formation of neutron stars with negligible natal kicks through electron-capture or accretion-induced supernova was also enabled. Another feature of these models is the inclusion of gravitational wave recoil kicks whenever two BHs merge (Baker et al. 2008; Morawski et al. 2018). Low birth spins were assumed for BHs, with values uniformly sampled between 0 and 0.1 (Fuller & Ma 2019). The orientation of the BH spin relative to the binary orbit was randomly distributed (Morawski et al. 2018). The full list of changes in the recent MOCCA code, together with their detailed description, can be found in Giersz et al. (2025).

Appendix B: BH subsystem evolution

The BHS plays a critical role in the development of 1P overconcentration and evolves differently from that in single-population GCs that are initially tidally underfilling. Figure B.1 shows the time evolution of the BHS-to-cluster mass ratio, highlighting five key phases.

thumbnail Fig. B.1.

Time evolution of the BH subsystem mass as a fraction of the total cluster mass. Labels 1-5 mark the distinct evolutionary phases described in Section B.

thumbnail Fig. B.2.

Time evolution of Lagrangian radii (10%, 50%, and 70%) for BHs, showing the progressive central concentration of the BHS

In phase 1, the BHS mass fraction rises rapidly due to BH formation from massive 1P stars and concurrent 1P mass loss as the cluster expands and overfills its tidal radius. In phase 2, the mass fraction drops following 2P formation, which significantly increases the total mass while contributing few BHs. The evolution beyond this point is shaped by BHS mass segregation and contraction (see Fig. B.2). In phase 3, this segregation fuels energy generation that drives enhanced stellar escape, causing the BHS mass fraction to rise again. The brief dip around 2 Gyr reflects the final stages of BHS formation and early BH ejections. Phase 4 marks a period of balanced evolution, with the cluster expanding and preferentially losing 1P stars. As the system is only mildly nTF, escape remains efficient, and the BHS mass fraction increases further. At 5.5 Gyr, the cluster migrates to a larger Galactocentric distance and becomes strongly nTF (phase 5), suppressing stellar escape and leading to gradual BHS depletion through internal interactions. While the present-day BHS mass fraction in our model (5-6%) is higher than values inferred for most Galactic GCs-which are typically below 1% (Weatherford et al. 2020; Dickson et al. 2024)-it is comparable to estimates for exceptional cases such as Omega Centauri (Zocchi et al. 2019; Baumgardt et al. 2019). For clusters like NGC 3201 and NGC 6101, the inferred BH mass fractions are typically lower, around 1-2%. However, we note that multi-mass modeling techniques used to derive these estimates can systematically underpredict BH mass fractions above 1%, particularly in dynamically young clusters (see Fig. 1 in Dickson et al. 2024). Although NGC 3201, similar to our model, has a large present-day half-mass radius and appears dynamically young, our simulation suggests that it underwent significant dynamical evolution prior to its migration to a larger Galactocentric distance.

Appendix C: Observational core radius

The observational core radius, as calculated in MOCCA, can vary substantially between output times, particularly at late stages of cluster evolution. This variability arises from the small number of bright evolved stars that dominate the central surface brightness, making the derived core radius sensitive to their spatial distribution. Figure C.1 shows the time evolution of the observational core radius in our model. While the long-term trend reflects the global dynamical state of the cluster, significant variations between output times are apparent –ranging from below 1 pc to over 6 pc. These fluctuations are not necessarily indicative of structural changes, but rather reflect changes in the central surface brightness profile caused by the positions of a few luminous stars. For example, at the output time shown in the main paper (13.002 Gyr), the observational core radius is 3.13 pc and the central surface brightness (CSB) is 610 L⊙/pc2. In the immediately following output (13.006 Gyr), these values shift to 2.07 pc and 817 L⊙/pc2, respectively-closer to the observational values reported for NGC 3201 in the Harris (1996, updated 2010) catalog, which lists a core radius of 1.85 pc and a CSB of 914 L⊙/pc2. Given these fluctuations, we emphasize that other structural parameters-such as total luminosity, cluster mass, and half-light radius-are more stable and provide a more reliable basis for comparison with observations. Additionally, the 3D core radius of the model at 13.002 Gyr is 1.22 pc, with a half-mass radius of 8.71 pc. At 13.006 Gyr, the 3D core radius is 1.77 pc, and the half-mass radius is 8.72 pc. These values are in excellent agreement with the parameters listed for NGC 3201 in the Baumgardt & Hilker (2018, updated 2024) catalog2.

thumbnail Fig. C.1.

Observational core radius as a function of time.

All Figures

thumbnail Fig. 1.

Cumulative number distributions of RGB stars for 1P (black line) and 2P (red line) as a function of distance scaled by Rhl for selected projected (2D) snapshots from 12 Gyr to 13.5 Gyr. 1P is more centrally concentrated than 2P only in the snapshot at 13 Gyr. The number distributions of both populations change significantly over time, and the overconcentration of 1P is a transient feature. The Rhl values, snapshot times, and the number of RGB stars are provided in the insets of each panel.
In the text
thumbnail Fig. 2.

Top left panel: Cumulative mass distributions of all stars for 1P (black line) and 2P (red line) as a function of distance scaled by Rhl in the projected (2D) snapshot at 13 Gyr. Top right panel: Cumulative number distributions of all stars for 1P (black line) and 2P (red line) as a function of distance scaled by Rhl in the snapshot at 13 Gyr. Bottom left panel: Cumulative number distributions of MS stars with masses greater than 0.4 M⊙ for 1P (black line) and 2P (red line) as a function of distance scaled by Rhl in the snapshot at 13 Gyr. Bottom right panel: Cumulative number distributions of RGB stars for 1P (black line) and 2P (red line) as a function of distance scaled by Rhl in the snapshot at 13 Gyr, shown up to a distance of Rhl.
In the text
thumbnail Fig. 3.

The lower panel: Escape time versus radial position for dynamically ejected present-day RGB progenitor stars with initial masses between 0.81 and 0.83 M⊙. Each point represents a present-day RGB progenitor that escaped due to a strong three-body (binary-single) or four-body (binary-binary) encounter. The x axis shows the escape time, and the y axis indicates the radial position within the cluster at which the escape-causing interaction occurred. Superimposed lines trace the evolution of the core radius (Rc, grey), half-mass radius (Rh, dark blue), half-light radius (Rhl, light blue), and tidal radius (Rt, brown). Violet points represent 1P escapers, and green points represent 2P escapers; filled symbols denote single stars, and open symbols denote binaries. Triangles indicate stars that escaped directly due to an encounter involving at least one BH, while squares indicate stars that had at least one such BH encounter during their dynamical history. Of the 159 1P escapers, 20 (12.6%) experienced BH interactions; among the 414 2P escapers, 96 (23.2%) were involved in BH interactions. The vertical dashed red line marks the 13 Gyr snapshot, where the 1P overconcentration is most evident. The top panel: Histograms of escape times for both populations, binned in 200 Myr intervals, highlighting the higher escape rate of 2P RGB progenitors from strong encounters.
In the text
thumbnail Fig. B.1.

Time evolution of the BH subsystem mass as a fraction of the total cluster mass. Labels 1-5 mark the distinct evolutionary phases described in Section B.
In the text
thumbnail Fig. B.2.

Time evolution of Lagrangian radii (10%, 50%, and 70%) for BHs, showing the progressive central concentration of the BHS
In the text
thumbnail Fig. C.1.

Observational core radius as a function of time.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.