© The Authors 2024

1 Introduction

Globular clusters (GCs) were previously believed to be simple stellar populations with a uniform age and chemical composition. However, thanks to extensive photometric and spectroscopic studies, it has become clear that these systems are more complex than previously thought and host multiple stellar populations (MSPs) characterized by differences in the chemical abundances of various light elements and, in some cases, also of iron (see, e.g., Gratton et al. 2019 for a review). This discovery has opened many new questions and challenged the traditional understanding of globular cluster formation and dynamical evolution.

The origin of MSP is still a matter of intense investigation and no consensus has been reached on the possible sources of processed gas necessary to explain the chemical composition of MSP and their star formation history. Some MSP formation models have proposed that all stars formed simultaneously, but some stars acquired an anomalous chemical composition by accreting processed gas produced by supermassive stars or massive binary stars (see, e.g., de Mink et al. 2009; Bastian et al. 2013, Gieles et al. 2018).

Other formation models suggest that after the formation of first-generation (FG) stars, a second episode of star formation takes place from gas enriched by the ejecta of FG polluters, and this results in the production of second-generation (SG) stars with different chemical abundances (e.g., enhancement in Na, Al, and N abundances and depletion in Mg, O, and C). In this scenario the origin of the gas that formed later generations of stars is of particular interest; possible candidates for the sources of polluted gas proposed in the literature include single and binary asymptotic giant branch (AGB) stars, rapidly rotating massive stars, massive interacting binaries, stellar mergers, and black hole accretion disks (see, e.g., Ventura et al. 2001; Decressin et al. 2007; de Mink et al. 2009; Elmegreen 2017; Krause et al. 2013; D’Ercole et al. 2008, 2010; Jěrábková et al. 2017; D’Antona et al. 2016; Breen 2018; Calura et al. 2019; Wang et al. 2020).

While the origin of MSP in globular clusters is not fully understood, different proposed scenarios based either on numerical hydrodynamical simulations (see, e.g., D’Ercole et al. 2008; Bekki 2010, 2011; Calura et al. 2019; Lacchin et al. 2022) or on general considerations (see, e.g., Gieles et al. 2018) typically share the prediction that SG stars tend to form more centrally concentrated compared to FG stars.

Although this difference in the initial spatial distributions of FG and SG stars is gradually erased during the cluster dynamical evolution, some clusters may still retain some memory of the structural properties imprinted by the formation process (see, e.g., Vesperini et al. 2013, 2021; Miholics et al. 2015; Sollima 2021 for some studies of the process of spatial mixing) and, indeed, several observational studies have found evidence of SG stars being more centrally concentrated than the FG population (see, e.g., Bellini et al. 2009; Lardo et al. 2011; Simioni et al. 2016; Dalessandro et al. 2019; Onorato et al. 2023; Lacchin et al. 2024). We also refer to Leitinger et al. (2023) for a study that found two clusters (NGC 3201, NGC 6101) where the FG is currently more centrally concentrated than the SG; however, Cadelano et al. (2024) performed an observational study of the structural and kinematic properties of one of the clusters, NGC 3201, providing support for scenarios in which the SG formed more centrally concentrated.

Differences between FG and SG stars are not limited to their spatial distributions but extend to their kinematic properties. Differences in the FG and SG kinematic properties may be imprinted during the formation process (see, e.g., Bekki 2010, 2011; Lacchin et al. 2022 for differences in the FG and SG rotation) or emerge during the cluster evolution (see, e.g., Vesperini et al. 2021 for differences between the anisotropy in the FG and SG velocity distribution). Evidence of kinematic differences have been found in several observational studies (see, e.g., Richer et al. 2013; Bellini et al. 2015, Cordero et al. 2017; Cordoni et al. 2020; Libralato et al. 2023).

The initial structural differences between the FG and SG populations may also have important implications for the survival and evolution of their binary stars. Hong et al. (2015, 2016, 2019) investigated the evolution of binary stars in MSP through direct NBODY simulations. As a consequence of the fact that the SG population is initially more centrally concentrated in a denser subsystem, SG binaries are either more easily disrupted compared to FG binaries or, for compact binaries, they evolve more rapidly toward more compact configurations. These investigations also predicted the presence of mixed binaries, composed of stars from different generations due to member exchange during strong interactions. Further detailed investigations of the dynamics of binary stars in multiple-population clusters and what they can reveal about the clusters’ initial structural properties have been carried out by Hypki et al. (2022) and Sollima et al. (2022).

On the observational side, the investigation of the binary populations in MSPs is still in its early stages. The findings of the few studies carried out to date are in general agreement with the predictions of numerical simulations concerning the fraction of FG and SG binaries and their variation with the clustercentric distance (D’Orazi et al. 2010; Lucatello et al. 2015; Dalessandro et al. 2018; Milone et al. 2020). Possible evidence of primordial differences in the fraction of FG and SG binaries have been suggested in the study by Kamann et al. (2020). The first evidence of the mixed binaries predicted by Hong et al. (2015) has been reported by Milone et al. (2020).

In this paper we continue our previous investigations of the dynamics of MSPs with an extensive survey of Monte Carlo simulations exploring a broad range of different initial conditions and shedding further light on the complex dynamics of MSP clusters. The goal of this paper is to carry out a general exploration of how the evolution of some of the fundamental MSP properties such as the fraction of SG stars and of some of the cluster’s structural properties depend on the cluster’s initial conditions. Although in some cases we compare our results to observations, we note that our analysis is not specifically aimed at reproducing the observed trends and distributions of the clusters’ observed properties; our general goal is rather to explore what initial conditions eventually lead to the present-day properties being within the range of those found in Galactic clusters and what their dynamical history is, and to expand the theoretical framework describing how various parameters affect the dynamics of MSP. The conclusions of this work will help to design the initial conditions for the next, more comprehensive MOCCA-SURVEY. A more extended survey and specific choices on the distribution of initial properties are required for a more detailed comparison with observations and will be the subject of future papers.

This paper is organized as follows. In Section 2 we describe the newest version of the MOCCA code, the initial conditions of the numerical simulations performed for this paper, and finally a short description of the data analysis software. Section 3 presents the results obtained with the MOCCA simulations through the context of the ratio of the number of objects from the SG to the total (N_SG/N_TOT). The MOCCA simulations are also briefly compared with the Milky Way (MW) GCs showing that MOCCA simulations are able to cover the observational ranges of cluster global properties, and thus are also good probes of the physical processes taking place between MSPs. In Section 4 we discuss the potential implications obtained from MOCCA simulations for the observational signatures of the multiple populations in GCs and some implications for the scenarios of their formation. Section 5 briefly summarizes the main findings of the paper.

2 Numerical simulations

This section presents the description of the MOCCA code, the simulations which were computed for this project, and the way the output data were analyzed.

2.1 MOCCA

This work is based on the numerical simulations performed with the MOCCA¹ Monte Carlo code (Giersz 1998; Hypki & Giersz 2013; Giersz et al. 2014; Hypki et al. 2022). MOCCA is a feature-rich advanced code that performs full stellar and dynamical evolution of real-size star clusters. Over the last few years it has been substantially updated and many new features have been added to the code, making it one of the most advanced and fastest codes in stellar dynamics able to simulate real-size star clusters up to Hubble time. The newest major additions include several features to support the study of the dynamics and stellar evolution of multiple stellar populations. A detailed description of the new MOCCA features is presented in Hypki et al. (2022).

MOCCA is able to follow the full dynamical evolution of MSPs and also the stellar evolution for different populations. Only mergers and mass transfers between stars from different populations are treated in a simplified way; the stars are marked as a mixed population because we do not provide procedures to accurately model the chemical mixing between two stars belonging to different populations.

The initial conditions explored in this paper and in our previous studies (Vesperini et al. 2021; Hypki et al. 2022) are inferred from the results of hydrodynamical simulations (see, e.g., D’Ercole et al. 2008; Calura et al. 2019) of multiple population formation from the ejecta of AGB stars and external pristine gas reaccreted by the cluster showing that SG stars form in a centrally concentrated subsystem embedded in a more extended FG system. As pointed out in Section 1, however, the prediction of such a spatial configuration is generally shared by other models based on different FG polluters.

Our simulations start with the FG and SG subsystems already in virial equilibrium and do not follow in detail the very early phases of SG formation. MOCCA allows stellar evolution to start for all populations at T = 0 or the stellar evolution of the SG can start after some time delay. The new physics that was added to the stellar evolution part of the code (additions to SSE/BSE code (Hurley et al. 2000, 2002), Belloni et al. (2017) is included. A comprehensive summary of the stellar evolution features of the MOCCA code can also be found in Kamlah et al. (2022). Strong dynamical interactions in MOCCA are performed with the FEWBODY code (Fregeau et al. 2004; Fregeau & Rasio 2007). The dissipative effects connected with tidal forces or gravitational wave radiation during dynamical FEWBODY interactions have not been taken into account yet; they are planned for the next version of the code.

2.2 Initial conditions

In this paper we extend the initial conditions explored in our previous studies and carry out a comprehensive investigation of the dependence of the evolution of the fraction of SG stars and the cluster structural parameters on the initial conditions. Following previous works (see Hong et al. 2016; Vesperini et al. 2021; Hypki et al. 2022, and references therein), the initial model contains a higher number of FG stars and the initial number ratio of the SG to FG stars is typically set between 0.33 to 0.38. It is also assumed that the SG is more centrally concentrated than the FG population, and we use the concentration parameter (conc_pop), defined as the ratio of the half-mass radii of SG to FG, to quantify the initial spatial differences between the two populations. Furthermore, in the initial model, the King (1966) concentration parameter (W₀) is specified separately for each population. For the SG, this was fixed to W_0,SG = 7 or 8, and for the FG population this parameter was varied. For FG stars, the initial zero age main sequence (ZAMS) masses were sampled between 0.08 and 150 M_⊙ . For SG stars, the upper limit for the ZAMS mass was set to 20 M_⊙ or 150 M_⊙. The ZAMS masses were sampled according to the Kroupa (2001) initial mass function (IMF).

We explore the evolution of models at a various fixed galac- tocentric distances, r_g , as indicated in Table 1. Similar to the models simulated in Hypki et al. (2022), the metallicity of both populations in all the simulated cluster models was set to Z = 0.001 (5% of Z_⊙). All models were simulated with an updated treatment for the evolution of massive stars (Tanikawa et al. 2020; Kamlah et al. 2022) 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 black holes (BHs) and neutron stars (NSs) 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 supernova was also enabled (Kamlah et al. 2022). Another feature of these models is the inclusion of gravitational wave (GW) recoil kicks whenever two BHs merge (Baker et al. 2008; Morawski et al. 2018). The magnitude of the GW recoil kick depends on the magnitude and orientation of the spins of BHs. In all the simulated models, low birth spins for BHs are assumed and are uniformly sampled values between 0 and 0.1 (Fuller & Ma 2019). The orientation of the BH spin with respect to the binary orbit is randomly distributed (Morawski et al. 2018).

The initial conditions are summarized in Table 1. Hereafter we refer to the models as the MOCCA-SURVEY. The table summarizes the initial conditions with all of the possible parameters. However, it is important to note that we did not compute all possible combinations of these parameters, but only a small subset of them. We computed over 200 MOCCA models for the purpose of this paper.

2.3 Data analysis

The data analysis for this paper was done via the BEANS² software (Hypki 2018). More precisely, the data analysis was performed in Apache Pig³, which is a high level language for the Apache Hadoop⁴ platform. BEANS allows users to have one script that is able to query various MOCCA simulations (from different surveys) in bulk and analyze them in a distributed way. Next, they can be shared and redistributed among the collaborators for further analysis. One example Apache Pig script is discussed in Hypki et al. (2022, Appendix B), and another is discussed in Appendix B of this paper.

3 Results

In this section we describe the results obtained from the MOCCA simulations while investigating some of structural properties of our models, and the evolution of N_SG/N_TOT. We note that our MOCCA models are not designed to reproduce the detailed properties of any specific globular cluster or the trends and correlations observed in the Galactic globular cluster system. In this paper our focus is on the study of the evolution of star clusters with multiple populations and how the fraction of SG stars and the clusters’ structural parameters depend on various initial properties. This will help us better constrain the initial conditions for a new MOCCA-SURVEY-2 containing a few thousand GC models. In some figures we plot the final values of some of the fundamental properties of our models along with the corresponding observed values of Galactic GCs; however, the goal is just to provide a general comparison of the range of final properties of our models with the corresponding observed values.

Fig. 1 Milky Way GCs coverage with MOCCA models for core (rc), half-light (rhl) radii, and the ratio between them (from top to bottom). The half-light and core radii for MW GCs are taken from Baumgardt et al. (2019). Only MOCCA models of at least 10 Gyr in age are shown.

3.1 Structural properties

In the three panels of Figure 1 we show the final values of core radius (r_c , distance at which surface brightness is equal to half of the central brightness), half-light radii (r_hl), and the ratio r_c/r_hl for all the MOCCA models that survived at least 10 Gyr as a function of the cluster’s final mass.

The three panels also show the observed values for Galactic GCs that were taken from Baumgardt et al. (2019). While our models were not specifically designed to replicate the properties of the Galactic GC system, it is interesting to note how tidally filling (TF) and tidally underfilling (TuF) MOCCA models both cover the range of values of the radii; in general they fall within the range of observed values. As already pointed out, exploration of a broader range of initial conditions is necessary to identify which models can evolve and reach values of the structural parameters not covered by our current survey. For example, larger values of rc and rhl may be attained by systems evolving at larger galactocentric distances than those studied in our survey and larger final values of the cluster masses may be obtained for systems initially more massive than those we have considered.

3.2 Evolution of N_SG/N_TOT ratios

In this section we explore the evolution of the N_SG/N_TOT ratio and the clusters’ half-light radii and how they depend on various initial properties of the cluster models. In Figure 2 (and Figure A.1 in Appendix A) we show the time evolution of these two quantities for various representative cases with initial conditions indicated in each panel of these figures. Overall the panels of Figures 2 and A.1 provide a view of the role of various parameters on the evolution of N_SG/N_TOT and show that in all cases the r_hl values are generally consistent with those found in Galactic clusters.

We divide the cases presented in these figures into two groups. The first group concerns TF, and the second TuF clusters. We discuss the role played by each of the parameters varied in our exploration.

The two panels in the Figure 2 show the evolution of the TF and TuF models with different total numbers of stars. The TF models all evolve very similarly reaching large values of N_SG/N_TOT of ~0.8 and half-light radii of ~4 pc consistent with those observed in many Galactic clusters. The most distinct difference is the longer dissolution times for models with larger N; this trend is simply the consequence of the higher masses and accordingly longer half-mass relaxation times and slower cluster evolution.

The TuF models also evolve toward final half-light radii of 3-4 pc, while the final N_SG/N_TOT are much lower than those found in the TF models and are just slightly lower than the initial values (~0.2). This is expected because these clusters undergo a weaker early loss of FG stars resulting in lower final values of N_SG/N_TOT, although they still fall in the range of values observed in Galactic clusters.

The MOCCA models for different galactocentric distances (r_g) are shown in Figure A.1 (first row). TF models with larger r_g have the same initial masses, but larger r_tid . All models rapidly lose FG stars during the early evolution dominated by the cluster’s expansion, which is driven by mass loss associated with stellar evolution; this results in a significant increase in N_SG/N_TOT approximately independent of the galactocen- tric distance. This result is consistent with the findings of Vesperini et al. (2021) who found that the final N_SG/N_TOT is mainly determined by the early evolutionary phases, while the effects of galactocentric distances only play the role of a much less important second parameter.

For our model with r_g = 2 kpc, the final N_SG/N_TOT is lower than for other models only because this cluster completely dissolves at the end of the simulation and the final stages of evolution are dominated by the mass loss from the center by dynamical interactions, preferably involving SG stars. For TuF models the N_SG/N_TOT does not significantly depend on r_g (from r_g = 2 kpc to 8 kpc). As expected, strongly underfilling models do not lose a significant fraction of FG stars, and the values of N_SG/N_TOT at the end of the simulation are similar to the initial values: they drop only slightly.

In the second row in Figure A.1 there are models with different binary fractions (fb). Models with different fb have the same total mass. TF models with higher fb have only slightly higher final N_SG/N_TOT. This is connected with faster mass segregation and larger energy generation in dynamical interactions by models with larger fb. Interestingly, it seems that fb = 0.1 (usually chosen in N-body simulations) results in a value for N_SG/N_TOT similar to that for larger fb. For TuF models fb does not have any visible influence on N_SG/N_TOT and only slightly on r_hl . This is connected with the fact that clusters behave as isolated models, and the structure of the cluster controls the efficiency of the central energy source.

The third row in Figure A.1 illustrates the results obtained by assuming different values for the maximum stellar mass (M_max) for the SG. The general expectation for single-population clusters is that a cluster with a higher M_max (or, more generally, with a larger initial fraction of massive stars) should undergo a stronger initial expansion because of larger amount of mass lost due to stellar evolution (see, e.g., Fukushige & Heggie 1995). However, the effect of varying M_max for the SG in a multiple-population cluster is more complex. Detailed investigation of the evolution of models with different M_max values reveals very interesting dynamical effects associated with the change in M_max. As expected for the TF model with higher SG M_max, the first phase of the cluster evolution connected with the stellar evolution leads to a larger mass loss and a more rapid increase in N_SG/N_TOT (the effect is more apparent for the case r_g = 4 pc and conc_pop =0.1, but is not shown in this paper). During the subsequent long-term evolution, however, the increase in N_SG/N_TOT slows down and even slightly decreases. The model with lower M_max , on the other hand, continues its evolution toward higher values of N_SG/N_TOT.

The different dynamical evolution of the two systems with different SG M_max can be explained by the underlying differences in the population of stellar mass BHs in the two systems. The model with SG M_max = 20 M_⊙ does not create a population of SG BHs, which are instead produced in the system with SG M_max = 150 M_⊙ . Thus, from the very early phases of its evolution the model with the higher SG M_max creates a dense core with a population of stellar mass BHs that act as an energy source. The model with SG M_max = 20 M_⊙, on the other hand, does not have very massive stars (and the BHs they form) in the centrally concentrated SG subsystem. BHs and massive stars have to segregate in the central regions from the more extended and less concentrated FG system. This process needs some time. It needs ∼2 Gyr to mass segregate 70% of BHs from FG. For the model with M_max = 150 M_⊙ only 1 Gyr is needed.

After mass segregation FG BHs will start to generate energy and support the cluster evolution. The longer timescale of mass segregation for the system with the smaller SG M_max leads to a stronger early expansion and to different spatial structures (Lagrangian radii for M_max = 20 M_⊙ are larger than for M_max = 150 M_⊙). The more extended structure of the systems with M_max = 20 M_⊙ leads to stronger mass loss of FG stars and higher N_SG/N_TOT and, more generally, to a more rapid cluster dissolution. The same process is much more profound and better visible for the case r_g = 4 pc and conc_pop = 0.1 (not shown in the paper), for which the model with M_max = 20 M_⊙ presents even larger increase in r_hl , and dissolves 4 Gyr faster than the model with M_max = 150 M_⊙.

A similar behavior and differences between models with different values of the SG M_max, although to a much smaller extent, are found for the TuF models. In the case of TuF models both models survive until the Hubble time.

The fourth row in Figure A.1 presents TF and TuF models for different values of the concentration parameter (conc_pop) that determines the relative spatial distributions of the FG and SG populations (see Section 2). For the TF models N_SG/N_TOT rapidly increases during the cluster early evolution for all the values of conc_pop; the value of N_SG/N_TOT at the end of the early evolution (at t ∼ 1–2 Gyr) is slightly higher for smaller values of conc_pop for which the higher concentration of the SG population leads to a stronger preferential loss of FG stars. The higher central densities of models with smaller values of conc_pop also affect the long-term evolution and the subsequent evolution of N_SG/N_TOT, but the general trend between N_SG/N_TOT and conc_pop is already set during the cluster’s early evolutionary phases.

For TuF models the loss of FG stars and the ensuing variation of N_SG/N_TOT are much milder and again the final values for N_SG/N_TOT are just slightly lower than the initial ones. Also for the models presented in these panels the final values of r_hl fall within the range of those found for Galactic globular clusters.

The fifth row in Figure A.1 shows the models for different W_0,FG which has the most profound influence on the N_SG/N_TOT ratio and dissolution time (t_diss). For TF models, the ratio of the FG half-mass to tidal radius increases for decreasing values of the W_0,FG values. This trend implies that GCs with smaller W_0,FG values lose FG stars more efficiently during the early cluster’s expansion, leading to a more significant increase in N_SG/N_TOT.

The process becomes less efficient for larger W_0,FG and for W_0,FG = 6 the value of N_SG/N_TOT increases only slightly in comparison to the initial values. This confirms earlier findings by Vesperini et al. (2021) that cluster models for which W_0,FG is too large (larger than 6) do not produce clusters with high present-day N_SG/N_TOT . We point out, however, that the inclusion of additional dynamical processes such as primordial gas expulsion and early tidal shocks might affect the structure of the clusters and lead to the efficient loss of FG stars, also for larger initial values of W_0,FG . For TuF models, the dependence of cluster parameters and N_SG/N_TOT on W_0,FG is practically negligible. Independently of W_0,FG , clusters have a lot of space to expand up to r_tid, and so the mass loss is very similar.

The last row in Figure A.1 shows models for different r_hFG (only TuF models). The r_tid radius is the same for all the models. These models nicely show that for larger r_hFG values and constant r_tid the ratio r_tid /r_hFG is increasing and coming closer to TF. Therefore, more FG stars escape, which leads to an increase in the ratio. Since all these models are tidally underfilling, GCs can undergo their initial expansion without immediately losing FG stars. and thus N_SG/N_TOT does not increase as much as TF models. Also in this case, the inclusion of the dynamical processes mentioned above (gas expulsion, early tidal shocks) may lead to a more significant loss of FG stars and an increase in the N_SG/N_TOT ratio.

For TuF models, we find an early slight decrease (about < 1%) in the N_SG/N_TOT ratio during the initial ~100 Myr. This is the result of the preferential ejection of SG stars and binaries due to dynamical interactions in the innermost regions (the main star loss mechanism in TuF clusters) where the SG is the dominant population.

Fig. 2 Evolution of NSG /NTOT ratios, together with half-light radii (rhl), for a set of MOCCA simulations. The figure consists of two series of MOCCA models: tidally filling (left, TF), and tidally underfilling (right, TuF). The main MOCCA initial parameters are summarized in the titles. In every plot there are 2–4 MOCCA models for which all initial parameters are the same, except one shown in the legend of each panel. In each panel, the left Y-axis shows the half-light radius (rhl (pc)), and the right Y-axis the NSG/NTOT ratio. The figure shows models differing in N, initial number of stars. The other varying parameters are: rg [kpc], galactocentric distance; fb, binary fraction; and Mmax [M⊙], maximum mass for the SG stars, and they are presented in an identical way in Figure A.1 in Appendix A.

3.3 The role of initial dynamical parameters on the evolution of clusters and their multiple stellar populations

As shown in the numerous past theoretical studies on the evolution of globular clusters, the dynamical history of these systems depends on various internal dynamical parameters and on the properties of the external tidal field of their host galaxies (see, e.g., Heggie & Hut 2003). The presence of MSPs with different dynamical properties significantly broadens the parameter space describing the possible initial properties of globular clusters and their subsequent evolution. The small survey of simulations carried out for this investigation allows us to start building a comprehensive picture of the possible dynamical paths followed by multiple-population clusters. Although a much larger number of simulations will be necessary to build a more complete picture, the simulations presented here provide a number of key indications of the role played by some of the parameters describing their initial structural properties.

Table 2 summarizes the results emerging from Figures 2 and A.1 and illustrates the role played by various parameters in determining the variation of N_SG/N_TOT, cluster mass at the Hubble time, or dissolution time. As expected t_diss increases as N (for TF, and TuF), r_g , W₀ , and M_max for TF models increase. In turn, in order to increase N_SG/N_TOT ratios one can increase fb for TF models, and r_hFG for TuF models.

4 Discussion

The main goal of this study is to provide an initial exploration of the role of various parameters in determining the evolution of a few key dynamical properties of globular clusters and those of their multiple stellar populations. This exploration is the first step toward the identification of the regions of the initial parameter space leading to properties generally consistent with the present-day observed properties of Galactic globular clusters. In the following sections we discuss further some of the results of the simulations introduced in this paper. In future studies we will extend our survey of simulations and carry out an investigation aimed at modeling the evolution of populations of globular clusters.

Fig. 3 Milky Way GCs coverage with MOCCA models for 2D NSG/NTOT (Milone et al. 2017a) ratios computed as a cumulative value for different radii (scaled to rhl). Every line correspond to one MOCCA simulation. The selected MOCCA models are the same as in Figures 2 and A.1. The NSG/NTOT ratios are colored according to whether they were initially TF or TuF. Only for a few models is there a large difference for NSG/NTOT ratios between small and large radial distances (for those closer to dissolution). For the majority of the models there is no significant difference for NSG/NTOT ratios between those measured in the very center and those measured around rhl radius. For details see Section 4.

4.1 Distribution of SG stars in the cluster

We start our discussion by focusing our attention on the radial variation of the fraction of SG stars. As discussed in the Section 1, several formation models and hydrodynamical simulations of SG formation predict that SG stars form more centrally concentrated in the inner regions and, during the subsequent early and long-term dynamical evolution, the strength of the initial spatial differences gradually decreases. Some clusters may retain some memory of these initial differences, while in others complete mixing may be reached. Observational studies often estimate N_SG/N_TOT in a limited portion of the cluster, and it is therefore important to establish how these estimates may be affected by radial variations in N_SG/N_TOT.

Figure 3 presents the N_SG/N_TOT ratio of MW GCs and the cumulative N_SG/N_TOT ratios computed for different radii in the MOCCA models. Each line corresponds to one MOCCA simulation at 12 Gyr (for models dissolving earlier that did not survive to this time, it is the profile calculated from the last saved snapshot that still holds at least 1% of the initial mass). The MOCCA models are those presented in Figures 2 and A.1. The cumulative N_SG/N_TOT ratios are computed for every MOCCA simulation for R_max from 0.1 r_hl , up to 10 r_hl . This figure clearly shows that despite the broad range of initial conditions explored and the fact that the models in our survey reach a variety of different degrees of spatial mixing, the values of N_SG/N_TOT measured within clustercentric distances similar to those usually covered by observational studies are representative of the global values. Moreover, as already discussed in the previous section, this figure further illustrates a clear dichotomy between initially TuF and TF models where the latter are generally characterized by larger final values of N_SG/N_TOT and that the values of N_SG/N_TOT found in our TF models falls within the range of those found in observations of Galactic GCs.

4.2 The range of values of N_SG/N_TOT

As shown in the various panels of Figures 2 and A.1, the values of N_SG/N_TOT in our models can reproduce the range of values observed in Galactic globular clusters (Milone et al. 2017b), while producing systems with half-light radii also generally consistent with those observed. Further investigation extending the survey of simulations to include a broader range of initial masses will allow us to study a populations of clusters and explore the trends between N_SG/N_TOT and other cluster properties.

We generally excluded models that formed an intermediatemass black hole (IMBH) from the analysis presented in this paper. However, some models do form them. We plan to study the formation and evolution of these IMBHs and their influence on the N_SG/N_TOT in a future work. In the newest MOCCA version, for example, we have IMBH seed BHs being formed as a result of a runaway merging scenario.

Our analysis shows that in order to reach values of N_SG /N_TOT consistent with those observed the FG populations must be initially TF or slightly TuF, and be characterized by initial values of W_0,FG smaller than about 5-6. More extreme initial central concentration (i.e., larger values of W₀) for the FG would lead to a modest early loss of FG stars and smaller increase in N_SG/N_TOT. Such values for W_0,FG are in agreement with those expected from the process of residual gas removal after FG is formed. Sudden gas removal leads to a much shallower central potential and a more extended FG characterized by small W₀ (Leveque et al. 2022).

For the requirement on the initial degree of tidal filling and the possible link between the initial and present-day structural properties, it is important to point out that, as shown in this study and previously in Vesperini et al. (2021), most of the evolution of N_SG/N_TOT occurs in the cluster early evolutionary phases with a much more modest increase in N_SG/N_TOT during the subsequent long-term evolutionary phase. This implies that the relevant dynamical requirements for the evolution of N_SG/N_TOT are to be considered in the context of the strength of the tidal field during the first phase of a cluster’s evolution. A number of studies (see, e.g., Meng & Gnedin 2022; see also El-Badry et al. 2016 for the possible effect of feedback-driven fluctuations in the gravitational potential of a galaxy in the early radial migration toward weaker tidal fields) have suggested that clusters experience a stronger tidal field (and stronger time variations and tidal shocks) in the first 1–2 Gyr of their evolution and later migrate to larger galactocentric distances and weaker tidal fields. The idea that short timescale cluster migration from gas-rich formation environments (e.g., via mechanisms frequent galaxy mergers) has been proposed as a mechanism for the long-term survival of GC progenitors (Kruijssen 2015; Forbes et al. 2018).

For the evolution of multiple-population clusters, in addition to consequences associated with the possible role of additional processes contributing to the early loss of FG stars (e.g., early tidal shocks), this migration implies that the required condition of tidally filling clusters would correspond to clusters with a more compact structure and smaller half-mass radii than those needed to be TF in weaker tidal fields. Regarding the evolution of the fraction of SG stars and the cluster’s mass, this could be a plausible pathway to support a modified scenario in which a sizable number of extended FG stars can be lost due to strong tidal stripping in the first 1–2 Gyr of cluster evolution in the initial stronger tidal field while the long-term evolution of the cluster’s mass driven by two-body relaxation would proceed at a slower rate for clusters migrating into a weaker tidal field (see Onorato et al. 2023 for a first investigation exploring the implications of such a transition for the evolution of the multiple populations in the massive Galactic cluster NGC 2419). The resulting masses, sizes, and N_SG/N_TOT in such a scenario incorporating a transition from a stronger to a weaker tidal field may naturally lead to properties consistent with those observed in present-day clusters. Further investigation of this scenario is currently in progress and will be presented elsewhere.

5 Conclusions and future work

In this paper we explored the evolution of multiple-population clusters for a broad range of initial conditions expanding those considered in previous studies. The exploration presented in this paper provides a more comprehensive picture of the evolution of a number of key properties of multiple-population clusters and will serve as the basis for future investigations and surveys aimed at building specific models for the properties of Galactic globular clusters and their multiple populations.

In our models we considered different initial numbers of stars, initial concentrations of the FG and SG populations, galac- tocentric distances, binary fractions, and upper mass limits on the initial mass function, as well as configurations in which the FG was tidally filling or tidally underfilling. This exploration allowed us to start shedding light on how a number of parameters depend on the initial conditions adopted. In particular, in this paper we started to explore how the cluster’s lifetime, total mass, and fraction of SG stars depend on the initial values assumed for those parameters and properties. The results are summarized in Table 2. In addition to the trends reported in Table 2, our conclusions can be summarized as follows:

• In agreement with previous studies, we find that in models starting with the FG tidally filling, N_SG/N_TOT can undergo a significant evolution reaching higher values falling in the range of those observed in Galactic globular clusters. Models with a FG initially tidally underfilling, on the other hand, do not lose a significant number of stars and retain values of N_SG/N_TOT similar to the initial values;

• In order for the clusters to undergo a significant increase in the N_SG/N_TOT ratio, the initial spatial distribution of the FG population modeled as that of a King model must have an initial value of the central dimensionless potential W₀ ∼ 5–6 or smaller;

• The N_SG/N_TOT ratio changes most noticeably during the first 1–2 Gyr of the cluster’s evolution and it does not change significantly during the subsequent evolution. The initial conditions and the environment in which a GC was born are thus likely to play crucial role in shaping the final values of N_SG/N_TOT ratios (see Figure 2);

• In most of the models we investigated for this paper, we find only mild differences (<0.1) between the value of N_SG/N_TOT calculated within the inner regions (e.g., within 0.1 r_hl) and the values calculated within 1–2 r_hl . In most cases values of N_SG/N_TOT calculated within 1-2 r_hl (a radial range typical of many observational studies) are representative of the global values for the entire cluster (see Figure 3);

• Many of the models and initial conditions explored in this paper produce final values of N_SG/N_TOT, masses, core and half-light radii overlapping with those observed in Galactic globular clusters (see Figures 1 and 2), and our survey has shed light on the range of initial conditions resulting in properties generally consistent with observations. We point out, however, that the goal of this paper was not to produce a complete model for the Galactic globular cluster system and the trends observed for the properties of the multiple populations. Additional simulations including clusters even more massive than those considered here are necessary for a more comprehensive investigation. Our future models will include additional ingredients and refinements; in particular, we will include the possible effects of a tidal field varying in time due to fluctuations in the cluster’s birth environment as well as a result the cluster’s migration from the site of formation to various galactocentric distances, the dynamical effects associated with a delay between the time of FG and SG formation, and tidal effects due to eccentric orbits. A significant extension of the survey presented here including these effects and exploring a broader range of initial conditions will be presented in future papers.

Data availability

All of the data used in this article can be shared on request. This includes the data from all MOCCA simulations like e.g. snapshots with stellar properties, positions, velocities of all stars and binaries, global properties of star cluster, and all derived data tables.

Acknowledgements

This research has been partially financed by the Polish National Science Centre (NCN) grant 2021/41/B/ST9/01191. EV acknowledges support from NSF grant AST-2009193. 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. 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. MOH acknowledges support by the Polish National Science Center grant 2019/32/C/ST9/00577. We thank the referee for all the comments and suggestions that helped us to improve the paper. The MOCCA code is open source⁵ for our collaborators. We are open to starting new projects, in which one could use already existing MOCCA simulations or start new ones. The BEANS⁶ software is open source and is freely available for anyone.

Appendix A Evolution of N_SG/N_TOT ratios

Fig. A.1 (cont. of the Figure 2) Every row, starting from the top, differs in: rg [kpc], galactocentric distance; fb, binary fraction; Mmax [M⊙ ], maximum mass for the SG stars, and concpop, concentration parameter. Every row, starting from the top, differs in: W0,FG, King parameter for FG; and rhFG [pc], half-mass radii for FG (thus, only TuF models).

Appendix B BEANS script

We present one of the scripts that was used while working on this paper. The whole data analysis was done with the BEANS software (see Section 2.3). We take this opportunity to show how one can analyze huge datasets (in our case astronomical data coming from numerical simulations) in an easy way using Apache Pig scripts. Other scientists might also be interested in using BEANS for their research.

A detailed description about specific Apache Pig keywords and instructions can be found in Hypki et al. (2022, Appendix B) or in the Apache Pig documentation.⁷ We advise reading it first. Here, only the main steps of the scripts are briefly described. The script computes cumulative profiles of FG and SG stars for all MOCCA simulations (used in this paper) for a number of selected timesteps for which there are available snapshot data. The comments in the Apache Pig scripts are the lines starting with the characters “- -” and they are used to describe code snippets below.

-- <1>  Reading data from snapshot files from MOCCA
--         simulations, from the Survey5, for a few
--         selected timesteps (0, 0.4, 1, 2, 5… Gyr)
snap = load ’datasets="mocca survey5"
                        tables="snapshot"
                        filter="(timenr == 0)
                        OR (tphys>400.0 tphys<430.0)
                        OR (tphys>1000.0 tphys<1030.0)
                        OR (tphys>2000.0 tphys<2030.0)
                        OR (tphys>5000.0 tphys<5030.0)
                        OR (tphys>7000.0 tphys<7030.0)
                        OR (tphys>10000.0 tphys<10030.0)
                        OR (tphys>12000.0 tphys<12030.0)"’
                        using BeansTable();
-- <2>  Getting from snapshots only some selected
--         columns: ID of the dataset (dsid), physical
--         time (tphys), radial distance (r),
--         population ID (pop), and we create bins
--         every r = 0.1 pc
snap = foreach snap generate
                        DSID(tbid) as dsid,
                        tphys,
                        r,
                        popId1 as pop,
                        FLOOR(r, 0.1) as bin;
-- <3>  Compute total number of FG and SG stars
--         for every MOCCA simulation, for every
--         timestep. First we group all data based
--         on these criteria (group … by), and
--         later we compute the numbers (foreach …
--         COUNT(…)) in every group. It will be
--         used to normalize the profiles later.
snapGr = group snap by (dsid, tphys, pop);
snapTotal = foreach snapGr generate
                        group.$0                        as dsid,
                        group.$1                        as tphys,
                        group.$2                        as pop,
                        COUNT(snap.pop)       as popCount;
-- <4>  Instructions in blocks <4>, <5> and <6>
--         concerns actual computation of the
--         profiles. They were divided into blocks
--         for clarity. Block <4> groups all data
--         by bins (with width r = 0.1 pc). Under
--         the variable (alias) ’snapGr’ all the data
--         are grouped by MOCCA simulations, physical
--         times, population id, and later by bins.
--         This is needed for the next instruction,
--         where one can compute how many FG and SG
--         are in every bin.
snapGr = group snap by (dsid, tphys, pop, bin);
snap = foreach snapGr generate
>                        group.$0 as dsid,
>                        group.$1 as tphys,
>                        group.$2 as pop, group.$3 as r,
>                        0COUNT(snap.r) as count;
-- <5>  Creating self-join between the results of
--         the previous code (alias ’snap’). In the
--         first step it is created a copy of alias
--         alias ’snap’. In the next step the actual
--         self-join operator is used. Then, in this
--         joined alias the rows are filtered and
--         the one comming from the ’snap2’ which have
--         larger radial distances are removed - this
--         is needed in the last <6> step to compute
--         cumulative values.
snap2 = foreach snap generate *;
snapJoined = join snap by (dsid, tphys, pop),
                        snap2 by (dsid, tphys, pop);
snapJoined = filter snapJoined
                        by snap::r <= snap2::r;
snapJoined = foreach snapJoined generate
                        snap::dsid                        as dsid,
                        snap::tphys                      as tphys,
                        snap::count                      as count,
                        snap2::r                           as r2,
                        snap2::pop                       as pop;
-- <6>  Computing cumulative profiles based on the
--         rows from the previous alias. The alias
--         is first grouped by the radial distance
--         ’r2’. Then, one can compute how many rows
--         there are inside such radius - these are
--         in fact cumulative values for both
-         populations FG, and SG.
snapJoinedGr = group snapJoined by (dsid, tphys,
                        pop, r2);
snap3 = foreach snapJoinedGr generate
                        group.$0 as dsid,
                        group.$1 as tphys,
                        group.$2 as pop,
                        group.$3 as r2,
                        SUM(snapJoined.count) as r2Count;
-- <7>  These two rows are simply joining the
--         results of the previous rows (with
--         cumulative values) with the rows which
--         store total number of FG and SG stars
--         (needed for normalization).
snap3Joined = join snap3 by (dsid, tphys, pop),
                        snapTotal by (dsid, tphys, pop);
snap3Joined = foreach snap3Joined generate
                        snap3::dsid                  as dsid,
                        snap3::tphy                 as tphys,
                        snap3::pop                   as pop,
                        snap3::r2                      as r,
                        snap3::r2Count            as count,
                        snapTotal::popCount  as popCount;
-- <8>  Storing the results in BEANS.
store snap3Joined into
                        ’name="snap cumulative survey5"’
                        using BeansTable();

References

All Tables

All Figures

	Fig. 1 Milky Way GCs coverage with MOCCA models for core (rc), half-light (rhl) radii, and the ratio between them (from top to bottom). The half-light and core radii for MW GCs are taken from Baumgardt et al. (2019). Only MOCCA models of at least 10 Gyr in age are shown.
In the text

Fig. 2 Evolution of NSG /NTOT ratios, together with half-light radii (rhl), for a set of MOCCA simulations. The figure consists of two series of MOCCA models: tidally filling (left, TF), and tidally underfilling (right, TuF). The main MOCCA initial parameters are summarized in the titles. In every plot there are 2–4 MOCCA models for which all initial parameters are the same, except one shown in the legend of each panel. In each panel, the left Y-axis shows the half-light radius (rhl (pc)), and the right Y-axis the NSG/NTOT ratio. The figure shows models differing in N, initial number of stars. The other varying parameters are: rg [kpc], galactocentric distance; fb, binary fraction; and Mmax [M⊙], maximum mass for the SG stars, and they are presented in an identical way in Figure A.1 in Appendix A.

In the text

Fig. 3 Milky Way GCs coverage with MOCCA models for 2D NSG/NTOT (Milone et al. 2017a) ratios computed as a cumulative value for different radii (scaled to rhl). Every line correspond to one MOCCA simulation. The selected MOCCA models are the same as in Figures 2 and A.1. The NSG/NTOT ratios are colored according to whether they were initially TF or TuF. Only for a few models is there a large difference for NSG/NTOT ratios between small and large radial distances (for those closer to dissolution). For the majority of the models there is no significant difference for NSG/NTOT ratios between those measured in the very center and those measured around rhl radius. For details see Section 4.

In the text

	Fig. A.1 (cont. of the Figure 2) Every row, starting from the top, differs in: rg [kpc], galactocentric distance; fb, binary fraction; Mmax [M⊙ ], maximum mass for the SG stars, and concpop, concentration parameter. Every row, starting from the top, differs in: W0,FG, King parameter for FG; and rhFG [pc], half-mass radii for FG (thus, only TuF models).
In the text