# Multiple stellar population mass loss in massive Galactic globular clusters

E. Lacchin, A. Mastrobuono-Battisti, F. Calura, C. Nipoti, A. P. Milone, M.

Meneghetti, E. Vanzella

► To cite this version:

E. Lacchin, A. Mastrobuono-Battisti, F. Calura, C. Nipoti, A. P. Milone, et al. Multiple stellar population mass loss in massive Galactic globular clusters. Astronomy & Astrophysics - A&A, 2024, 681, 10.1051/0004-6361/202347268. insu-04821684

HAL Id: insu-04821684

https://insu.hal.science/insu-04821684v1

Submitted on 6 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Multiple stellar population mass loss in massive Galactic globular clusters

E. Lacchin<sup>1,2,3,4,5</sup>, A. Mastrobuono-Battisti<sup>3,6,7</sup>, F. Calura<sup>1</sup>, C. Nipoti<sup>2</sup>, A. P. Milone<sup>4,8</sup>, M. Meneghetti<sup>1</sup>, and E. Vanzella<sup>1</sup>

Received 23 June 2023 / Accepted 25 September 2023

## ABSTRACT

The degree of mass loss, that is the fraction of stars lost by globular clusters, and specifically by their different populations, is still poorly understood. Many scenarios of the formation of multiple stellar populations, especially the ones involving self-enrichment, assume that the first generation (FG) was more massive at birth than now in order to reproduce the current mass of the second generation (SG). This assumption implies that, during their long-term evolution, clusters lose around 90% of the FG. We tested whether such strong mass loss could take place in a massive globular cluster orbiting the Milky Way at 4 kpc from the centre that is composed of two generations. We performed a series of *N*-body simulations for 12 Gyr to probe the parameter space of internal cluster properties. We derive that, for an extended FG and a low-mass SG, the cluster loses almost 98% of its initial FG mass and the cluster mass can be as much as 20 times lower after a Hubble time. Furthermore, under these conditions, the derived fraction of SG stars, *f*<sub>enriched</sub>, falls in the range occupied by observed clusters of similar mass (~0.6–0.8). In general, the parameters that affect the highest degree of mass loss are the presence or absence of primordial segregation, the depth of the central potential, *W*<sub>0,FG</sub>, the initial mass of the SG, *M*<sub>SG</sub><sup>ini</sup>, and the initial half-mass radius of the SG, *r*<sub>h,SG</sub>. Higher *M*<sub>SG</sub><sup>ini</sup> have not been found to imply higher final *f*<sub>enriched</sub> due to the deeper cluster potential well which slows down mass loss.

**Key words.** methods: numerical – globular clusters: general – stars: kinematics and dynamics – Galaxy: evolution – Galaxy: kinematics and dynamics – Galaxy: disk

## 1. Introduction

In the last decades, an increasing number of observations have revealed the presence of multiple stellar populations (MPs) within globular clusters (GCs), a discovery that has revolutionised our view of these stellar systems (Gratton et al. 2019). Stars belonging to distinct populations differ in their light element abundances (such as C, N, O, Na, Mg, and Al), while they share – at least in the bulk of GCs – the same iron content. These variations are well defined and linked by anticorrelations such as the C–N, Na–O, and Mg–Al ones (Piotto et al. 2005; Carretta et al. 2009; Milone et al. 2017; Gratton et al. 2019; Masseron et al. 2019; Marino et al. 2019). In particular, within the same GC we can distinguish between stars sharing the same chemical composition of the field ones (O-rich and Na-poor), labelled as first population, and O-poor and Na-rich stars classified as the second population. However, different populations not only differ in their chemical abundances, but also in their structural and kinematical properties, suggesting a deep connection between the origin of the chemical imprints in MPs and the formation and the subsequent dynamical evolution of the whole cluster.

Although the long-term dynamical evolution these systems have undergone is gradually erasing the structural and kinematical differences that MPs had at birth, dynamically young clusters may retain some memory of the original differences between MPs until the present day, in particular in their outskirts (Vesperini et al. 2013). By means of observational data analysis, supported by *N*-body models, Dalessandro et al. (2019) showed the tight connection between the relative degree of concentration of different populations and the evolutionary stage of the cluster. In particular, larger differences in the spatial radial distributions between distinct populations are found in clusters that are dynamically young and have experienced lower mass loss. The difference between the two populations is not restricted to the spatial distribution; observational studies have revealed that second-population stars are, in some clusters, characterised by a more radially anisotropic velocity distribution (Richer et al. 2013; Bellini et al. 2015, 2018; Libralato et al. 2019, 2023), a more rapid rotation (Lee 2015, 2017, 2018; Cordero et al. 2017; Dalessandro et al. 2019; Kamann et al. 2020; Cordini et al. 2020; Szigeti et al. 2021), a lower fraction of binaries (D’Orazi et al. 2010; Lucatello et al. 2015; Milone et al. 2020), and a more central concentration

<footnote n="1">¹ INAF – OAS, Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, Via Gobetti 93/3, 40129 Bologna, Italy
e-mail: elena.lacchin@inaf.it</footnote>

<footnote n="2">² Dipartimento di Fisica e Astronomia, Università di Bologna, Via Gobetti 93/3, 40129 Bologna, Italy</footnote>

<footnote n="3">³ GEPI, Observatoire de Paris, PSL Research University, CNRS, Place Jules Janssen, 92195 Meudon, France</footnote>

<footnote n="4">⁴ Dipartimento di Fisica e Astronomia “Galileo Galilei”, Univerisità di Padova, Vicolo dell’Osservatorio 3, 35122 Padova, Italy</footnote>

<footnote n="5">⁵ INFN – Padova, Via Marzolo 8, 35131 Padova, Italy</footnote>

<footnote n="6">⁶ Department of Astronomy and Theoretical Physics, Lund Observatory, Box 43, 221 00 Lund, Sweden</footnote>

<footnote n="7">⁷ Max Planck Institute for Astronomy, Königstuhl 17, 69117 Heidelberg, Germany</footnote>

<footnote n="8">⁸ Istituto Nazionale di Astrofisica – Osservatorio Astronomico di Padova, Vicolo dell’Osservatorio 5, 35122 Padova, Italy</footnote>
(Norris & Freeman 1979; Sollima et al. 2007; Lardo et al. 2011; Milone et al. 2012; Richer et al. 2013; Cordero et al. 2014; Simioni et al. 2016; Dondoglio et al. 2021) than the first population.

In addition to the different degrees of variations in the struc-
tural and kinematical properties between MPs, more massive
clusters are generally found to host a higher fraction of second-
population stars (up to 90%) than lower mass ones (down to
30–40%, and even 10% in the Magellanic Clouds; see e.g.,
Milone & Marino 2022; Milone et al. 2017; Zennaro et al. 2019;
Dondoglio et al. 2021). This quantity is the result of a complex
combination of formation history and evolutionary effects, since,
due to the differences between the two populations at birth, they
will experience distinct dynamical evolution and therefore dis-
tinct mass-loss rates, which implies a change of the fraction of
the second population with time.

Despite a large amount of observational and theoretical stud-
ies providing new insights into the chemical and kinematical
properties of MPs, a clear understanding of how globular clus-
ters were formed has still not been reached (Renzini et al. 2015;
Bastian & Lardo 2018; Gratton et al. 2019). One of the crucial
points deals with the origin of the processed material and the
consequent formation of stars with an ‘anomalous’ chemical
composition. Many scenarios have been suggested in order to
tackle this issue proposing different sources of the processed
gas such as asymptotic giant branch (AGB) stars (D’Ercole et al.
2008, 2016; Bekki et al. 2017; Calura et al. 2019), fast-rotating
massive stars (Decressin et al. 2007), massive stars (Elmegreen
2017), supermassive stars (Denissenkov & Hartwick 2014;
Gieles et al. 2018), massive interacting binaries (de Mink et al.
2009; Bastian et al. 2013; Renzini et al. 2022), black hole accre-
tion discs (Breen 2018), and stellar mergers (Wang et al. 2020).
Nevertheless, so far, none of these scenarios has been able
to reproduce all the available observational constraints, and
therefore further and more thorough developments are required
(Renzini et al. 2015; Bastian & Lardo 2018).

The physical processes modulating the mass loss in stel-
lar clusters are manifold. Firstly, two-body relaxation was
found to gradually set up a Maxwellian velocity distri-
bution, which leads loosely bound stars to overcome the
cluster escape velocity (Ambartsumian 1938; Spitzer 1940).
Later, Chernoff & Weinberg (1990) followed the evolution
of multi-mass clusters with a tidal cut-off, driven by two-
body relaxation and stellar evolution mass loss. They found
that the combination of these two processes leads to a
stronger mass loss than the sum of the two independent
contributions. The dynamical evolution of a stellar clus-
ter is, however, affected by many other factors, such as
binarity (Tanikawa & Fukushige 2009; Fujii & Portegies Zwart
2011), tidal fields (Baumgardt & Makino 2003), gravitational
and tidal shocks (Gnedin & Ostriker 1997; Vesperini & Heggie
1997), mass segregation (Baumgardt et al. 2008; Vesperini et al.
2009; Haghi et al. 2014), and the presence of dark remnants
(Contenta et al. 2015; Banerjee & Kroupa 2011; Giersz et al.
2019). All these quantities are, however, known for several
present-day clusters, but not for star-forming clusters, posing
challenges in setting the initial conditions for simulated clus-
ters. Such uncertainty on the initial values also affects many
other parameters, such as the initial mass of the cluster (and
its radial distribution), which would be vital to understand-
ing how clusters form and dynamically evolve. Indirect deriva-
tions can be obtained starting from clusters’ present-day mass
and fraction of enriched stars. Assuming that different pop-
ulations are also distinct generations, and therefore that GCs

have undergone self-enrichment, it is possible to define a first
generation (FG) composed of normal stars and a second gen-
eration (SG) whose stars possess the peculiar chemical compo-
sition. If the mass of the FG is assumed to be comparable to
the present-day mass of GCs, one ends up with a mass released
by the FG polluters that is much lower than the mass of SG
stars observed today, which leads to the ‘mass budget prob-
lem’. To overcome this problem, it is generally assumed that
the cluster, and therefore the FG, was much more massive –
between 5 and 20 times – (Decressin et al. 2007; D’Ercole et al.
2008; Schaerer & Charbonnel 2011; Cabrera-Ziri et al. 2015) at
its birth; however, during the evolution, most of the FG stars (up
to ~95%) were lost, so the observed relative number of SG and
FG is still reproduced.

Few attempts have been carried out to determine whether,
during the long-term evolution, clusters are able to lose a
significant fraction of FG stars and then reproduce, after a
Hubble time, the observed clusters' features. Pioneering work
on the topic was done by D'Ercole et al. (2008), who per-
formed a series of N-body simulations in the AGB frame-
work, concluding that a cluster with a more concentrated SG
generation loses a substantial number of FG stars, at vari-
ance with the SG ones, deriving fractions of main sequence
(MS) stars $f_{\text{MS}} = N_{\text{SG,MS}}/N_{\text{FG,MS}}$ in agreement with obser-
vations. Later, Bastian & Lardo (2015) showed that no match
was found when combining the observational data with the
results of the N-body studies of Baumgardt & Makino (2003)
and Khalaj & Baumgardt (2015), concluding that neither gas
expulsion nor the effect of tidal fields could lead to the present-
day fraction of SG.

By means of Monte Carlo simulations, Vesperini et al. (2021) and Sollima (2021) studied the dynamical evolution of a cluster composed of two populations and a mass of ~$10^6 M_\odot$. Similarly to D’Ercole et al. (2008), they found that the cluster loses more FG stars and reaches, after 13 Gyr, the typical values of SG fraction observed in present-day GCs. Similar results were also obtained by Sollima et al. (2022), which focused on the binary fractions of the populations. They concluded that the present-day SG binary fraction can be used to constrain the initial concentration of SG stars, providing a relation between the initial size of the SG and total cluster mass.

From E-MOSAICS cosmological simulations, Reina-Campos
et al. (2018) explored the impact of dynamical cluster disruption
of multiple stellar populations deriving the degree of mass loss
and the fraction of enriched stars as a function of cluster mass,
galactocentric distance, and metallicity. They found discrepan-
cies with observations and therefore concluded that mass loss is
unlikely to have a strong impact on shaping the present-day GCs.
They also derived that to reconcile the observations a signifi-
cantly larger half-mass radius has to be assumed at birth, and a
higher initial SG fraction than the currently adopted ones would
be necessary.

Although the fraction of enriched stars is a very strong con-
straint widely used to compare simulated clusters with observed
ones, other fundamental pieces of information can be extracted
from the unbound stars (Arunima et al. 2023). Larsen et al.
(2012) found that around 1/5 of the metal-poor stars in the
Fornax dwarf spheroidal galaxy belong to the four GCs, mean-
ing that these GCs could have been, at most, five times more
massive at their birth, posing a strong upper limit on the frac-
tion of stars that could have been lost by GCs. Besides Fornax,
GCs have been found to resemble the Galactic ones (Larsen et al.
2014), and therefore they could have shared a common origin
and evolution; stars initially belonging to Fornax GCs could
have been lost in the early phases, therefore loosening the con-
straints about the degree of mass loss suffered by the Fornax
GCs (Khalaj & Baumgardt 2016). Based on the stellar chemical
composition, several studies have been carried out aimed
at determining the contribution that GCs could have given
to the formation of the Galactic halo (Carretta et al. 2010;
Martell & Grebel 2010; Martell et al. 2011, 2016; Ramírez et al.
2012). Recently, Koch et al. (2019) have analysed the spectra
of halo field giant stars. They found that 2% of the stars in
the sample show the ‘anomalous’ chemical composition typ-
ical of SG stars, in agreement with the previous investiga-
tions. In addition, they derived that 11% of the stars in the
Galactic halo were formed in GCs. This quantity is, however,
strongly affected by the adopted mass-loss rate in the early
phases and the number of completely dissolved clusters, reach-
ing up to 40–50% when assuming a mass loss factor, that is
the ratio between initial and final cluster mass, greater than ten
(Vesperini et al. 2010). Both the fraction of field SG stars and
that of field GC stars are extremely precious, providing fur-
ther constraints to the models, not only on cluster scales, but
also on larger ones, to understand how the Galaxy assembly
proceeded.

In this paper, we aim to derive the degree of mass loss in the
two different stellar components in order to determine whether
there are combinations of initial parameter values that can lead
to a significant mass loss – as with the one required to solve
the mass budget problem – and that spawn final clusters com-
patible with the observed GCs. We performed a series of direct
N-body simulations to follow the long-term evolution of a glob-
ular cluster with an initial mass of M ~ 10⁷ M⊙ and com-
posed of two populations taking into account stellar evolution,
the tidal effects of the Galactic potential, and primordial segrega-
tion. Although GC mass loss has been explored in the past, only a
few works have modelled a GC composed of more than one pop-
ulation (i.e., Vesperini et al. 2021; Sollima 2021; Sollima et al.
2022, with Monte Carlo codes). Our simulations are among the
first of this kind performed with a direct N-body code in the
literature, together with the ones of D’Ercole et al. (2008) and
Hénault-Brunet et al. (2015), where, however, a lower mass clus-
ter has been considered. The cluster is composed of two stellar
populations, and it is assumed to orbit the Milky Way (MW).
From the results derived by Calura et al. (2019), also confirmed
by Lacchin et al. (2022), the fraction between FG and SG is
larger than the present-day ones assumed to solve the mass bud-
get problem. Here, we aim to quantify the mass-loss factors that
are needed to reproduce the observed SG fraction for such a mas-
sive GC. The cluster is located in the disc of the Milky Way
and, therefore, is intended to represent a cluster belonging to
the in situ population. There are various reasons why studying
disc GCs is important. First, the distribution of metal-rich MW
GCs are more concentrated and flatter than the metal-poor com-
ponent, and generally, they are associated with the thick disc
and bulge populations (Armandroff & Zinn 1988; Armandroff
1989; Zinn 1985; Minniti 1995; Côté 1999; Van Den Bergh
2003; Bica et al. 2006, 2016). In addition, disc GCs, which
now constitute almost one-third of the total MW GCs (Harris
2010), could have been much more numerous in the past. Field
stars showing GC-like features have been discovered in the
inner Galaxy (Schiavon et al. 2017; Fernández-Trincado et al.
2022), a detection supported by simulations showing that
tidal effects in the inner regions of MW-like galaxies could
have gradually destroyed disc GCs, decreasing their popula-
tion (Renaud et al. 2017). Lastly, kinematic heating due to sev-
eral accretion events could have also deprived the disc of GCs,

which would now be part of the inner halo (Kruijssen 2015;
Di Matteo et al. 2020).

The paper is organized as follows. In Sect. 2, we describe
the model we adopted and the novelty introduced in the present
work. Section 3 deals with the results we obtained for our sets
of simulations. In Sect. 4, we discuss the outcomes of the sim-
ulations and compare them with the literature and observations.
Finally, we draw our conclusions in Sect. 5.

## 2. Models and method

We studied the internal evolution of a series of cluster mod-
els, including either one or two stellar populations, firstly with-
out star formation and stellar evolution and then including these
ingredients to study their effect on mass loss. In the next sub-
sections, we provide details on the initial setups adopted for
our simulations, illustrating the assumptions we adopted depend-
ing on the characteristics of the model and the phenomena that
we explored.

2.1. *Description of the code*

We ran our simulations using an updated version of NBSymple
(Capuzzo-Dolcetta et al. 2011)<fnref n="1" />, a direct and symplectic N-
body code parallelised on GPUs. Several versions of this
code are available and have been used to study various
aspects of GCs evolution in the Galactic potential (see e.g.,
Mastrobuono-Battisti et al. 2012, 2019; Sollima et al. 2012;
Leigh et al. 2014). This new version includes a star formation
and stellar evolution routine for the first and second generations
separately. The softening length adopted to avoid close encoun-
ters is equal to the mean interparticle distance within the half-
mass radius as done by D’Ercole et al. (2008) and it is calculated
separately for the first and second generations. Our choices for
the softening length and the number of particles used to represent
the clusters are motivated by the computational limitations, as
we aimed at 12 Gyr long simulations. We have tested the depen-
dence of our results on the choice of the softening length by
performing a simulation with a softening length ten times lower
than the mean interparticle distance, with the time step mod-
ified accordingly, as in Mastrobuono-Battisti et al. (2012). We
obtain that the cluster properties are only very weakly affected,
in particular the fraction of SG stars, f_enriched. Our simulations
are not intended to model close encounters and hard binaries; to
account for these effects, we would need a substantially higher
number of particles, smaller softening length, and shorter time
step, which would make 12 Gyr simulations computationally
unfeasible. However, we recall that our simulations start when
first-generation massive stars have already exploded. In addi-
tion, in order to avoid iron pollution, second-generation stars
are assumed to be composed of only low and intermediate-mass
stars. Therefore, massive interacting binaries are not included
in our models either because they have already evolved or they
never form. Lower mass binaries, neglected in our model but
expected to be present in real systems, would increase the num-
ber of ejected stars (Küpper et al. 2008). However, while binaries
heat up the system, a smaller softening makes the cluster more
compact, limiting its ability to lose stars. In any case, all these
collisional effects are expected to influence the mass loss much
less than those related to stellar evolution, which are accounted
for in our models. The mass lost due to stellar evolution is

<footnote n="1">¹ See https://github.com/alessandramb/NBSymple for the basic version of the code.</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Illustration of the adopted Galactic and GC models. The left panel provides a schematic representation of the Galaxy model adopted for the simulations, including both the thick and thin disc plus the halo, while the bulge is modelled as part of the Galactic disc. The right panel illustrates the GC model. The simulated cluster is located at 4 kpc from the Galactic centre and is assumed to orbit it on the plane of the disc. Credit: Alessandra Mastrobuono-Battisti.</figcaption></figure>

instantaneously removed, and, therefore, energy and angular momentum are not conserved.

## 2.2. Galactic potential model

Our Galactic model consists of a dark matter halo with both a thin and thick disc, as shown in the left panel of Fig. 1 (see also Mastrobuono-Battisti et al. 2019). The functional forms of these components are taken from Allen & Santillan (1991), with the parameters from Model II of Pouliasis et al. (2017), which is aimed at reproducing the actual MW. Such a model is able to reproduce various observables, such as the rotation curve, thin and thick disc scale lengths and heights, and stellar density in the solar neighbourhood. The bulge is considered as part of the Galactic disc, so it is not represented as an independent component (Di Matteo 2016). The mass assumed for the halo is $2.07 \times 10^{11} M_{\odot}$, with a scale height of 14 kpc. The thick disc has a mass of $3.91 \times 10^{10} M_{\odot}$, with radial and vertical scale length of 2 kpc and 800 pc, while the thin disc has a mass of $3.68 \times 10^{10} M_{\odot}$, a radial scale length of 4.8 kpc, and a scale height of 250 pc. With the adopted analytic Galactic model, we do not take into account the dynamical friction, which would be self-consistently included if the MW had been modelled as an N-body system, thus composed of stellar and dark-matter particles. To assess the importance of this neglected process, we estimated, through Eq. (8.13) of Binney & Tremaine (2008), that for a satellite of mass $8.6 \times 10^6 M_{\odot}$ the timescale of dynamical friction at 4 kpc from the Galactic centre is $t_{df} = 80$ Gyr/ln$\Lambda$, where ln $\Lambda$ is the Coulomb logarithm, assuming $M(<4$ kpc) = $4 \times 10^{10} M_{\odot}$ for the mass of the MW within 4 kpc derived through the model of Pouliasis et al. (2017). Despite the spherical approximation<fnref n="2" />, this result suggests that dynamical friction is negligible for the systems we are modelling.

## 2.3. Single stellar population models

Our single population clusters are modelled using the King (1966) profile, adopting different values for the half-mass radius, $r_h$, and for the dimensionless central potential, $W_0$ (Binney & Tremaine 2008), which varies between two and

seven. In this way, we explore the behaviour of both loose and dense clusters. The cluster initial mass is, for all models, $10^7 M_{\odot}$, with an initial metallicity of $Z = 0.001$, since, as explained before, we aim to model very massive clusters, such as the ones in Calura et al. (2019) and Lacchin et al. (2021, 2022). The initial metallicity of the cluster is $Z = 0.001$, as in Calura et al. (2019). The initial positions and velocities of the particles, in the absence of the external gravitational galactic potential, are derived using the software NEMO, through the *mkking* routine (Teuben et al. 1995). Our clusters are single-mass models (i.e., all stellar particles have the same mass) and are represented either with $N = 102400$ or $N = 25600$ particles, meaning that each particle is significantly more massive than a star. This choice is due to current computational limitations in running direct N-body simulations of systems with $10^6$ or more particles for a time span of 12 Gyr.

The adoption of single-mass stars has implications on the mass loss suffered by the clusters. We did not model mass segregation from binaries and massive stars remnants, which would sink into the innermost regions of the cluster, eventually favouring the loss of low-mass stars. However, the major driver of mass loss is stellar evolution, which also induces mass loss due to the shallowing of the potential, so we would not expect that the main conclusions of this work to be significantly different if one considered simulations accounting for the presence of a mass spectrum.

Each cluster in our simulation starts with an actual mass that is slightly lower than $10^7 M_{\odot}$. Since the simulation starting point is set after the explosion of FG core-collapse supernovae (i.e., at a time equal to $t_0 = 30$ Myr), we remove 16% of the initial cluster mass reaching $M_0 = 8.4 \times 10^6 M_{\odot}$, which mimics the effect of the death of massive stars. This value is the mass return fraction due to the evolution of stars with a mass larger than $8 M_{\odot}$ for the Kroupa (2001) initial mass function (IMF) with a final-to-initial mass relation from Agrawal et al. (2020).

Due to this mass removal, the system goes out of virial equilibrium and will expand to return to an equilibrium state. In our best model, this relaxation leads to a half-mass radius increase of 7% and no significant change in the mass loss, since the relaxation is taking place in the inner region, while the outskirts are very weakly affected.

We explored both clusters with and without primordial mass segregation. Primordial mass segregation has been found to have

<footnote n="2">² See Bonetti et al. (2021) for dynamical friction calculations for disc structures.</footnote>
<table><caption>**Table 1.** Models with neither stellar evolution nor SG stars run in this work.</caption><thead><tr><th>Model <sup>(a)</sup></th><th>W<sub>0</sub></th><th>r<sub>h</sub><br>(pc)</th><th>Segregation</th><th>Number of particles</th><th>Softening length<br>(pc)</th><th>f<sub>mass loss</sub></th></tr></thead><tbody><tr><td>n7N</td><td>7</td><td>23</td><td>N</td><td>25 600</td><td>1.65</td><td>0.44</td></tr><tr><td>n7Y</td><td>7</td><td>23</td><td>Y</td><td>25 600</td><td>1.65</td><td>0.58</td></tr><tr><td>n5N</td><td>5</td><td>37</td><td>N</td><td>25 600</td><td>2.57</td><td>0.50</td></tr><tr><td>n5Y</td><td>5</td><td>37</td><td>N</td><td>25 600</td><td>2.57</td><td>0.57</td></tr><tr><td>n2N</td><td>2</td><td>60</td><td>N</td><td>25 600</td><td>4.20</td><td>0.72</td></tr><tr><td>n2Y</td><td>2</td><td>60</td><td>Y</td><td>25 600</td><td>4.20</td><td>0.79</td></tr><tr><td>N2N</td><td>2</td><td>60</td><td>N</td><td>102 400</td><td>2.65</td><td>0.57</td></tr></tbody></table>

**Notes.** The initial mass of each simulated cluster is $M_0 = 8.4 \times 10^6 M_\odot$, after the removal of 16% of its mass. <sup>(a)</sup>Model name: n or N = small or large number of particles + W<sub>0</sub> + N or Y = non-segregated or segregated. Columns: 1) name of the model; 2) adimensional central potential parameter W<sub>0</sub> of the FG; 3) half-mass radius of the FG; 4) primordial segregation of the FG (N = non-segregated, Y = segregated); 5) number of particles N<sub>tot</sub>; 6) softening length; 7) mass-loss fraction defined as $f_{\text{mass loss}} = (M_{\text{ini}} - M_{\text{fin}})/M_{\text{ini}}$.

<table><caption>**Table 2.** Models with simplified stellar evolution and without SG stars run in this work.</caption><thead><tr><th>Model <sup>(a)</sup></th><th>W<sub>0</sub></th><th>r<sub>h</sub><br>(pc)</th><th>Segregation</th><th>Number of particles</th><th>Softening length<br>(pc)</th><th>f<sub>mass loss</sub></th></tr></thead><tbody><tr><td>n7Ne</td><td>7</td><td>23</td><td>N</td><td>25 600</td><td>1.65</td><td>0.63</td></tr><tr><td>n7Ye</td><td>7</td><td>23</td><td>Y</td><td>25 600</td><td>1.65</td><td>0.74</td></tr><tr><td>n5Ne</td><td>5</td><td>37</td><td>N</td><td>25 600</td><td>2.57</td><td>0.70</td></tr><tr><td>n5Ye</td><td>5</td><td>37</td><td>Y</td><td>25 600</td><td>2.57</td><td>0.79</td></tr><tr><td>n2Ne</td><td>2</td><td>60</td><td>N</td><td>25 600</td><td>4.20</td><td>1.00</td></tr><tr><td>n2Ye</td><td>2</td><td>60</td><td>Y</td><td>25 600</td><td>4.20</td><td>1.00</td></tr><tr><td>N2Ne</td><td>2</td><td>60</td><td>N</td><td>102 400</td><td>2.65</td><td>1.00</td></tr></tbody></table>

**Notes.** The initial mass of each simulated cluster is $M_0 = 8.4 \times 10^6 M_\odot$, after the removal of 16% of its mass. <sup>(a)</sup>Model name: n or N = small or large number of particles + W<sub>0</sub> + N or Y = non-segregated or segregated + e = with stellar evolution. Columns: same as in Table 1.

a significant effect on the cluster mass loss due to the cluster expansion in response to the massive star mass loss, happening preferentially at the cluster centre (Vesperini et al. 2009; Haghi et al. 2014). In case the clusters are mass segregated, we use the software McLuster (Küpper et al. 2011) to calculate the radius comprising all the massive stars in a primordially segregated model. We then remove the mass that is lost due to the explosion of stars more massive than $8 M_\odot$ within this radius, keeping a King profile for the density. All the models orbit the Galaxy in the plane of the disc at a galactocentric distance of 4 kpc, as shown in the right panel of Fig. 1. This is the same distance assumed by D’Ercole et al. (2008), albeit they used a simpler model, including a point-like mass located at the galaxy centre. It is worth mentioning that in our model, the mass enclosed within a radius of 4 kpc, derived integrating the density distribution, is $4.0 \times 10^{10} M_\odot$. This is in very good agreement with the mass assumed for the point-mass galaxy potential of D’Ercole et al. (2008). The clusters are tidally filling, i.e. their tidal radius, $r_t$, is equal to the distance at which the cluster potential and the Galactic potential have the same value (von Hoerner 1957; Baumgardt & Makino 2003; Webb et al. 2013). As the tidal radius is fixed to 200 pc, the core radius and half-mass radius of each of the models vary depending on the value of the W<sub>0</sub> parameter.

### 2.3.1. Without long-term stellar evolution

We initially modelled clusters hosting a single stellar population, which are only affected by dynamical effects (i.e., not considering any long-term stellar evolution effect). To start with the same

cluster mass, we removed the 16% of the initial mass, as we did in all the other models. The details on the models are reported in Table 1.

### 2.3.2. Adding long-term stellar evolution

In our second set of models, we still have only one stellar population, but we considered the effects of long-term stellar evolution. This is done through a mass return fraction taken from the relation between remnant mass and progenitor of Agrawal et al. (2020, case METISSE with MESA of Fig. 7) given by

$$m_{\text{loss}}(t) = m_{\text{p}}(t=0)[b_0 + b_1 \log(t) + b_2 \log^2(t)], \quad (1)$$

where $m_p(t=0)$ is the mass of the particle before the removal of the 16% of the mass due to massive star winds and SN explosions, $b_0 = 0.329420$, $b_1 = 0.0379353$, $b_2 = -0.002760463$, and t is expressed in Gyr. At t = 0.03 Gyr, the mass lost is the 16% of the whole mass, which is the mass we removed to take into account the death of massive stars. It is worth noting that at metallicity $Z > 10^{-3}$, as in our case, the variations in the fractional cumulative mass loss (expressed by the “returned fraction”) are of the order of a few percent (Vincenzo et al. 2016). The parameters adopted for these models can be found in Table 2.

### 2.4. Two stellar population models

In our third set of simulations, we finally added the SG, embedded inside the FG component. All the particles in the system have the same mass, for a total number of particles N<sub>tot</sub> = 102 400.
Both components are spherical and represented by King (1966) models. The FG component, modelled with $N_{\text{FG}}$ particles, is a King model with $W_{0,\text{FG}}$ ranging between two and seven with a total initial mass of $M_{\text{FG}}^{\text{ini}} = 10^7 M_{\odot}$ and a tidal radius of 200 pc, to mimic a tidally filling system. As before, we built the initial positions and velocities using the software NEMO, through the `mkking` routine (Teuben et al. 1995).

The SG component is a King model where we vary $W_{0,\text{SG}}$ from five to seven, its mass from $M_{\text{SG}}^{\text{ini}} = 7 \times 10^5 M_{\odot}$ to $3 \times 10^6 M_{\odot}$ (as a consequence, the number of SG particles $N_{\text{SG}} = N_{\text{tot}} - N_{\text{FG}}$ will also change from 7877 to 26 947, respectively), and the half-mass radius from 1 to 6 pc. We varied the SG mass to test different initial SG fractions. We vary the velocity dispersion of the SG as well, to explore the effect of this parameter on the cluster's mass-loss rate. We ran models with different values of the central velocity dispersion, equal to 0, $10 \text{ km s}^{-1}$, and to the velocity dispersion of the generated King model, where for the first two values we rescaled the velocities derived for the King model. In the third case, the SG is in equilibrium as an isolated system, while in the other cases it is radially out of equilibrium and tends to collapse and readjust after a phase of violent relaxation.

As before, the mass of the FG at the beginning of the simulation is slightly lower than its initial mass, since 16% of the mass is removed due to the explosions of core-collapse supernovae. After that, the FG starts to evolve dynamically as FG stars lose their mass due to stellar evolution, with a cumulative mass return fraction given by Eq. (1).

The SG appears after 10 Myr from the beginning of the simulation (i.e., at a time of $t_1 = 40 \text{ Myr}$) and grows its mass at a constant star formation rate of $0.05 M_{\odot} \text{ yr}^{-1}$ (for $M_{\text{SG}}^{\text{ini}} = 3 \times 10^6 M_{\odot}$; $0.01 M_{\odot} \text{ yr}^{-1}$ for $M_{\text{SG}}^{\text{ini}} = 7 \times 10^5 M_{\odot}$) for a total of 60 Myr (see Calura et al. 2019). To avoid the contribution of SG massive stars, which would chemically pollute the AGB ejecta with iron, for example, we assume – as is generally done in the AGB scenario – a truncated SG IMF composed only of stars with masses smaller than $8 M_{\odot}$ (see D’Ercole et al. 2010; Bekki 2019). SG stars are kept fixed with respect to the cluster centre of density while they are forming. Once the total initial mass of the SG is reached, they start to evolve dynamically. After an additional 30 Myr, the SG has accumulated enough mass and its stars start evolving following a cumulative mass return fraction law of the same shape of Eq. (1), but rescaled by a factor of 1.27133. The time evaluation of $m_p$ is shifted as well, due to the later formation of the SG, and corresponds to the time at which the SG stopped growing in mass.

The mass is added or removed in equal measure from each star particle in the relevant component of the cluster.

# 3. Results

In this section, we present the results obtained from our simulations. First, we describe the outcomes of models where only the FG is modelled and stellar evolution is not taken into account. Secondly, we report the results of the simulations assuming stellar evolution but still with the FG component only. Lastly, the outcomes of the simulations with both stellar evolution, FG, and SG components are described.

## 3.1. Models with a single stellar population

### 3.1.1. Models without stellar evolution

We first studied the long-term evolution of a massive cluster, $M_{\text{ini}} = 10^7 M_{\odot}$, composed only of FG stars. At the beginning

of the simulation, the stellar mass is equal to $M_0 = 8.4 \times 10^6 M_{\odot}$, which represents the mass of low- and intermediate-mass stars plus the remnants of the massive ones left in the system after massive stars have exploded.

In Table 1, we summarise the main parameters of our models together with the resulting mass-loss fraction, $f_{\text{mass loss}} = (M_{\text{ini}} - M_{\text{fin}})/M_{\text{ini}}$, with $M_{\text{fin}}$ being the final mass of the whole cluster at the end of each simulation. Stellar evolution is not taken into account for the moment. Figure 2 shows the evolution of the cluster mass, in the top left, and of the normalised mass loss, in the bottom left. As one can expect, the shallower the potential well, and therefore the lower the values of $W_0$, the greater the mass loss at the end of the simulation. However, model n7Y, which is characterised by $W_0 = 7$, is losing more mass than models with lower $W_0$ values in the first few million years. The higher concentration coupled with initial segregation is responsible for this behaviour, which also affects other models described throughout the section. Initial segregation generally leads to a stronger mass loss since segregated systems will have, after the death of massive stars, a larger $r_h$, making the system less bound. Apart from varying physical parameters, we also changed the number of particles from 25 600 to 102 400; hence, when more particles are used, as in N2N, the cluster loses less mass as a result of the longer relaxation time.

### 3.1.2. Models with stellar evolution

Figure 2 shows, on the right, the mass (top) and normalised mass-loss (bottom) evolution for the models listed in Table 2, where stellar evolution is taken into account. For comparison, the initial conditions we have adopted here are the same as for the models without stellar evolution described above.

As expected, mass removal due to stellar evolution leads to a shallower cluster potential well and, therefore, spurs subsequent mass loss in the form of lost stars. For models with $W_0 = 2$, the addition of stellar evolution leads to the dissolution of the cluster after ~3 Gyr due to the initially shallow potential. In all other cases, stellar evolution is less catastrophic, even though the final mass of the cluster is significantly smaller, from one-third to half, with respect to the case without stellar evolution. As before, initially segregated clusters suffer a stronger mass loss than non-segregated ones.

## 3.2. Models with second-generation stars

In Table 3, we report the initial conditions and the final values of masses, half-mass radii, fraction of enriched stars belonging to the final cluster (e.g., considering only bound stars), the fraction of unbound SG stars, and the central density for the thirteen simulations we performed, taking into account the stellar evolution and including the SG. It has to be stressed that, as for the previous models, the reported initial value for $r_{h,\text{FG}}$ is not the half-mass radius at the time of FG formation, but after the gas expulsion and violent relaxation phases, when the system is considerably more extended than at its formation. During these phases, the half-mass radius of a cluster can increase by a factor of three or four (Lada et al. 1984; Baumgardt & Kroupa 2007). Its exact value depends on many parameters, e.g., the IMF, the star formation efficiency, the gas and stellar density, and the gas expulsion timescale. Such large initial radii are also confirmed by observations of star-forming clusters at high redshift, where systems extending for several tens of parsecs have been detected. Further discussion regarding the scale radius of star-forming stellar clusters is reported in Sect. 4.3.
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Mass (top) and mass loss (bottom) as a function of time. The left column represents the models with the FG only and without stellar evolution, listed in Table 1, while in the right column the models with only FG but where stellar evolution is taken into account, listed in Table 2, are shown.</figcaption></figure>

In all models, we assumed that the SG is initially more centrally concentrated than the FG, and, at the end of the simulations and at different degrees all clusters still show this configuration. We varied several parameters in order to determine their effects on the evolution of the system focusing on the fraction of SG stars, which has been determined for several GCs, both in the Milky Way and in external galaxies (Milone et al. 2020; Dondoglio et al. 2021). Here we calculate it as the mass of bound SG particles over the mass of all the bound particles, $f_{\text{enriched}}(<r) = M_{\text{SG}}(<r)/M_{\text{tot}}(<r)$, both within the half-mass radius and for the entire cluster. In Sect. 3.2.1, we discuss the caveats that need to be considered when comparing the theoretical value to the observed one.

In all models, clusters start with the same FG mass, and given the possible non-self-similarity of the SG formation (e.g., different initial SG fraction for clusters of different initial mass, see Yaghoobi et al. 2022), we did not study the trend of $f_{\text{enriched}}$ with cluster mass here; this will be addressed in a future work.

First, we studied models with different values of the central velocity dispersion, equal to 0 km s⁻¹, 10 km s⁻¹ and the velocity dispersion of the generated King model. We find that such a quantity weakly affects the evolution, and the resulting clusters possess very similar fractions of SG stars. For this reason, all the subsequent simulations were performed assuming the same velocity dispersion distribution, corresponding to the King one.

As found for the previous models, with and without stellar evolution and without the SG, a segregated cluster loses a larger amount of mass than a non-segregated one. Comparing models 5N7M3.7K and 5Y7M3.7K, we derive that the FG loses significantly more mass in the segregated system, while the SG mass is weakly affected, since we only imposed the segregation on the FG. As a consequence, the $f_{\text{enriched}}$ is higher for the segregated system, especially within the half-mass radius.

We then varied the $W_{0,SG}$ parameter, but we found very weak effects on the evolution of the system, both in terms of mass loss by the two populations, $f_{\text{enriched}}$, half-mass radii, and central density (see the pairs 5Y7M1K – 5Y5M1K, 2N7m1K – 2N5m1K, 2N7m6K – 2N5m6K).

The parameters whose variation has a stronger and more complex impact on the mass loss and $f_{\text{enriched}}$ are, instead, $W_{0,FG}$, $M_{\text{SG}}^{\text{ini}}$, and the initial $r_{h,SG}$. In general, analogously to what has been derived for the simulations without SG, the greater the concentration of the FG, and therefore the larger the $W_{0,FG}$ value, the less mass lost by the whole cluster, as well as by the two populations separately, at the end of the simulation. While for $W_{0,FG} = 5$, a low mass of the SG equal to $7 \times 10^5 M_{\odot}$ (model 5Y7m1K) leads to a mild FG mass loss and a small final $f_{\text{enriched}}$ (0.16), for $W_{0,FG} = 2$, the same initial SG mass (e.g., model 2N7m1K) allows a loss of more than 95% of the FG mass. Consequently, when $W_{0,FG} = 2$, the final FG mass is one order of
<table><caption>**Table 3.** Models with both simplified stellar evolution and SG stars.</caption><thead><tr><th rowspan="2">Model <sup>(a)</sup></th><th colspan="6">Initial</th><th colspan="6">Final</th></tr><tr><th>W<sub>0,FG</sub><br>(pc)</th><th>r<sub>h,FG</sub><br>(pc)</th><th>Segr</th><th>M<sub>SG</sub><br>(10<sup>5</sup> M<sub>☉</sub>)</th><th>W<sub>0,SG</sub><br>(pc)</th><th>r<sub>h,SG</sub><br>(pc)</th><th>M<sub>FG</sub><br>(10<sup>5</sup> M<sub>☉</sub>)</th><th>M<sub>SG</sub><br>(10<sup>5</sup> M<sub>☉</sub>)</th><th>f<sup>r<sub>h,tot</sub></sup><sub>enriched</sub> (f<sub>enriched</sub>)</th><th>f<sup>lost</sup><sub>SG</sub><br>(pc)</th><th>r<sub>h,FG</sub><br>(pc)</th><th>r<sub>h,SG</sub><br>(pc)</th><th>r<sub>h,tot</sub><br>(pc)</th><th>log(ρ<sub>c</sub>)<br>(M<sub>☉</sub> pc<sup>-3</sup>)</th></tr></thead><tbody><tr><td>5N7M3.70</td><td>5</td><td>37</td><td>N</td><td>30</td><td>7</td><td>3.7</td><td>35.4</td><td>18.0</td><td>0.44 (0.34)</td><td>0.16</td><td>25.0</td><td>12.5</td><td>19.9</td><td>3.61</td></tr><tr><td>5N7M3.710</td><td>5</td><td>37</td><td>N</td><td>30</td><td>7</td><td>3.7</td><td>33.6</td><td>18.5</td><td>0.49 (0.35)</td><td>0.15</td><td>25.5</td><td>10.1</td><td>18.9</td><td>4.08</td></tr><tr><td>5N7M3.7K</td><td>5</td><td>37</td><td>N</td><td>30</td><td>7</td><td>3.7</td><td>36.3</td><td>19.8</td><td>0.49 (0.35)</td><td>0.14</td><td>24.1</td><td>9.34</td><td>17.3</td><td>4.11</td></tr><tr><td>5Y7M3.7K</td><td>5</td><td>37</td><td>Y</td><td>30</td><td>7</td><td>3.7</td><td>26.3</td><td>19.3</td><td>0.58 (0.42)</td><td>0.13</td><td>27.4</td><td>9.22</td><td>18.2</td><td>4.41</td></tr><tr><td>5Y7M1K</td><td>5</td><td>37</td><td>Y</td><td>30</td><td>7</td><td>1</td><td>22.9</td><td>15.7</td><td>0.46 (0.40)</td><td>0.16</td><td>27.8</td><td>22.4</td><td>25.4</td><td>2.36</td></tr><tr><td>5Y5M1K</td><td>5</td><td>37</td><td>Y</td><td>30</td><td>5</td><td>1</td><td>21.0</td><td>15.5</td><td>0.48 (0.42)</td><td>0.15</td><td>27.7</td><td>22.4</td><td>25.2</td><td>2.36</td></tr><tr><td>5Y7m1K</td><td>5</td><td>37</td><td>Y</td><td>7</td><td>7</td><td>1</td><td>21.9</td><td>4.42</td><td>0.23 (0.16)</td><td>0.03</td><td>32.1</td><td>20.2</td><td>29.9</td><td>2.15</td></tr><tr><td>2N7M1K</td><td>2</td><td>60</td><td>N</td><td>30</td><td>7</td><td>1</td><td>12.6</td><td>14.1</td><td>0.58 (0.52)</td><td>0.15</td><td>24.4</td><td>19.9</td><td>21.8</td><td>2.34</td></tr><tr><td>4N7M4.5K</td><td>4</td><td>45</td><td>N</td><td>30</td><td>7</td><td>4.5</td><td>27.8</td><td>18.8</td><td>0.58 (0.40)</td><td>0.13</td><td>29.3</td><td>7.84</td><td>19.4</td><td>4.90</td></tr><tr><td>2N7m6K</td><td>2</td><td>60</td><td>N</td><td>7</td><td>7</td><td>6</td><td>3.02</td><td>3.33</td><td>0.56 (0.52)</td><td>0.04</td><td>13.2</td><td>10.8</td><td>11.8</td><td>2.55</td></tr><tr><td>2N5m6K</td><td>2</td><td>60</td><td>N</td><td>7</td><td>5</td><td>6</td><td>2.93</td><td>3.50</td><td>0.59 (0.54)</td><td>0.03</td><td>12.9</td><td>10.2</td><td>11.3</td><td>2.45</td></tr><tr><td>2N7m1K</td><td>2</td><td>60</td><td>N</td><td>7</td><td>7</td><td>1</td><td>1.89</td><td>2.00</td><td>0.59 (0.52)</td><td>0.05</td><td>10.9</td><td>6.19</td><td>8.32</td><td>4.98</td></tr><tr><td>2N5m1K</td><td>2</td><td>60</td><td>N</td><td>7</td><td>5</td><td>1</td><td>2.15</td><td>2.09</td><td>0.57 (0.49)</td><td>0.05</td><td>10.7</td><td>6.11</td><td>8.25</td><td>4.47</td></tr></tbody></table>

**Notes.** The initial FG mass of each simulated cluster is $M_0 = 8.4 \times 10^6 M_☉$, after the removal of 16% of its mass. <sup>(a)</sup>Model name: $W_{0,FG} + Y$ or $N$ = with or without segregation + $W_{0,SG} + m$ pr M = low or high initial $M_{SG} + \text{initial } r_{h,SG} + \sigma_{0,SG}$ (K = King velocity dispersion). Columns: 1) name of the model; 2) $W_0$ of the FG; 3) half-mass radius of the FG; 4) primordial segregation of the FG (N = non-segregated, Y = segregated), 5) mass of the SG; 6) $W_0$ of the SG; 7) half-mass radius of the SG; 8) total mass of FG bound stars; 9) total mass of SG bound stars; 10) fraction of SG stars within the half-mass radius of the whole cluster (fraction of SG stars of the whole cluster); 11) fraction of SG stars among unbound stars $f_{\text{lost,SG}} = M_{\text{lost,SG}}/M_{\text{lost,tot}}$, where $M_{\text{lost,SG}}$ and $M_{\text{lost,tot}}$ are the mass of unbound SG and the total mass of unbound stars, respectively; 12) half-mass radius of the bound FG; 13) half-mass radius of the bound SG; 14) half-mass radius of the whole cluster (i.e., only bound stars); 15) central density.

magnitude lower than in the model with $W_{0,FG} = 5$, and the cluster reaches a final $f_{\text{enriched}}$ of 0.59 slightly lower than the typical values (~0.6–0.8; see Milone et al. 2017; Dondoglio et al. 2021) observed in GCs of the same mass. Similarly, model 2N7M1K loses more FG mass than 5Y7M1K, reaching a final SG fraction of 0.58.

Interestingly, however, for fixed $W_{0,FG} = 2$, variations of $M_{\text{SG}}^{\text{ini}}$ lead to significantly different final masses and $r_{h,tot}$, but $f_{\text{enriched}}$ slightly changes (models 2N7M1K – 2N7m1K). In model 2N7M1K, a higher SG mass implies a larger initial SG fraction; so, even though the cluster has lost almost an order of magnitude less mass, and consequently has a cluster radius more than double, the final SG fraction is very similar to the one of model 2N7m1K. On the other hand, for fixed $W_{0,FG} = 5$, variations of $M_{\text{SG}}^{\text{ini}}$ lead to similar FG masses but significantly different SG ones (see 5Y7M1K – 5Y5m1K). Therefore, the $f_{\text{enriched}}$ values differ by almost a factor of 2. These differences in the evolution suggest that there is not a positive correlation between the initial and final values of $f_{\text{enriched}}$. This is also visible in Fig. 3, where models assuming the same $W_{0,FG}$ do not always follow the same evolution. Therefore, clusters with an initially higher SG fraction do not necessarily have a higher final one.

A similar behaviour can be found when changing the initial $r_{h,SG}$. While for clusters with a massive initial SG a smaller $r_{h,SG}$ leads to a stronger SG expansion and lower $f_{\text{enriched}}$ (see models 5Y7M3.7K – 5Y7M1K), for an initial low-mass SG, clusters with smaller initial $r_{h,SG}$ are more compact at the end of the simulations, with final SG fractions similar to the ones of models starting with larger $r_{h,SG}$ (see models 2N7m6K – 2N7m1K and 2N5m6K – 2N5m1K). Therefore, an initially more compact SG does not imply a lower SG mass loss.

Comparing all the models with observed GCs more quantitatively, with a particular focus on the final masses and half-mass radii, it is clearly visible that models assuming $M_{\text{SG}}^{\text{ini}} = 3 \times 10^6 M_☉$ and $W_{0,FG} = 5$ produce a cluster with a mass almost two

times greater than $\omega$ Centauri, the most massive globular cluster known to date (Baumgardt & Hilker 2018). The very massive SG prevents a significant mass loss in these systems, and therefore the SG fraction is also too small (0.45–0.50) in comparison with the observations. The only exception, in terms of SG fraction, is model 5Y7M3.7K, where the combination between FG segregation and large initial $r_{h,SG}$ implies a greater loss of FG while poorly affecting the SG. Similar results are obtained decreasing $W_{0,FG}$, like in models 2N7M1K and 4N7M4.5K, where a stronger FG mass loss is taking place, leading to a higher final $f_{\text{enriched}}$. Moreover, all models with a massive SG are characterised by a fraction of lost SG stars, with respect to the total unbound mass, $f_{\text{SG}}^{\text{lost}}$, of ~0.15. This is between three and five times larger than what is found for a low-mass SG.

On the other hand, clusters with $M_{\text{SG}}^{\text{ini}} = 7 \times 10^5 M_☉$ lose much more mass, which depends on the $W_{0,FG}$. As highlighted above, for low $W_{0,FG}$, clusters undergo a strong mass loss, especially in the FG, resulting in final clusters with $M \sim 4–6 \times 10^5 M_☉$, in agreement with present-day ones.

In Fig. 3, we show the evolutions of the models 2N5m1K, 2N5m6K, 2N7M1K, 5Y7m1K, 5Y7M3.7K, and 5Y7M1K to highlight the effects of changing $W_{0,FG}$, $M_{\text{SG}}^{\text{ini}}$ and $r_{h,SG}$. Models 2N5m1K and 2N5m6K are the ones undergoing the strongest mass loss, and consequently their $r_h$ and $f_{\text{enriched}}$ also suffer deep changes during the evolution in the opposite direction; while $r_{h,tot}$ decreases by about six times after 12 Gyr, $f_{\text{enriched}}$ increases by almost the same amount.

Interestingly, comparing the results obtained with SG and $W_{0,FG} = 2$ with the ones with the same $W_{0,FG}$ but without SG, we clearly see that the formation of a concentrated SG within a shallow FG prevents the disruption of the clusters. An SG, even if not very massive, located at the centre of the system, is enough to strengthen the potential well of the cluster, decreasing the potential energy of the particles that will become more bound.
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Evolution of half-mass radius (top left), mass (top right), SG mass fraction of the whole cluster (bottom left), and mass loss (bottom right) for some of the models in Table 3 reported in the legend. The thick lines in the top panels represent the whole cluster and the thin SG.</figcaption></figure>

### 3.2.1. Model 2N7m1K

We focus here on the analysis of the model 2N7m1K, which is the one whose final fraction of SG, defined as $f_{\text{enriched}}(r) = M_{\text{SG}}(<r)/M_{\text{tot}}(<r)$, is in better agreement with current observations that find fractions between 50% and 80% (Milone et al. 2017). In Fig. 4, we show the distribution of the two populations on the left and the projected surface density on the right, both at 1 and 12 Gyr, with a zoomed-in image centred in the centre of mass of the system in the small inset in the middle of each panel. Already at 1 Gyr, the interaction of the initially spherical cluster with the Milky Way tidal field causes a distortion of the system, which develops two significant tidal tails – one leading (internal) and one trailing (external) – departing from the centre of the disrupting cluster. Most of the stars in the tails belong to the FG, which, after 12 Gyr, is still more extended than the SG, as it was at the beginning of the simulation. Later, the tails lengthen, reaching the main body, and two concentric circles appear in the maps. At 12 Gyr, the tails are dominated by the FG, while only 5% of the particles belong to the SG. Due to the intense mass loss suffered by the FG, the cluster is significantly more compact, as is clearly shown in the surface density maps. It is also more dominated by SG stars. Even though the system has been

highly distorted, its central region preserves a spherical shape after 12 Gyr.

As expected, the two populations, which were spatially and kinematically different at the beginning, move towards a mixing. This is spatially highlighted by the change in the density profiles and, in turn, in the half-mass radii. Such a behaviour can be clearly seen in Fig. 5, where we display the density profiles for the two populations together with the one of the whole cluster, compared with the profiles at the end of the simulation. While the SG is always dominant in the centre and its central density does not vary significantly over time, the FG undergoes a notable change in its profile. The FG increases its central density and decreases its half-mass radius; this results from the loss of stars in the outskirts due to the interaction with the Galactic tidal field. Overall, the central density of the cluster is in good agreement with the ones derived in present-day GCs (Baumgardt & Hilker 2018), while its half-mass radius of 8.3 pc is slightly larger than the ones of GCs with mass ~4 × 10⁵ M☉ (McLaughlin & van der Marel 2005; Krumholz et al. 2019). It loses almost 98% of FG stars, as generally predicted to match the observed SG fraction, with a final mass-loss factor of about 20, which is significantly smaller than the one reported in various other studies.
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Two-dimensional maps at 1 and 12 Gyr for the model 2N7m1K. Left panel: distribution of FG (cyan) and SG (magenta). Right panel: mass surface density of whole cluster. At the centre of both panels, we show a zoomed-in image focused on the centre of mass of the cluster.</figcaption></figure>

Although the two populations become more and more mixed with time, at the end of the simulation they are still well distinguished, with the SG being more concentrated than the FG. Such a difference in the shape of the profiles of the two populations affects the radial fraction of SG stars, $M_{SG}/M_{tot}$, which decreases with the distance from the Galactic centre, as shown in Fig. 6. This means that the fraction of SG stars within a fixed radius $r$, $f_{enriched}(r)$, is not flat all over the cluster, but it decreases from around 0.7 in the centre, down to 0.54 when considering the whole cluster. Such a decrease has been observed in Milky Way GCs such as 47 Tucanae and NGC 5927 (Milone et al. 2012; Cordero et al. 2014; Dondoglio et al. 2021; Jang et al. 2022), and in simulations (D’Ercole et al. 2008). Regarding $f_{enriched}$, it is important to stress that the observational values of this quantity are rarely calculated for the whole cluster, but, due to the

limited field of view (e.g., Milone et al. 2017), they often refer to the fraction of enriched stars within the inner regions of a cluster (typically between the centre and $r \sim 0.5–1.5r_h$). It is therefore important to compute $f_{enriched}$ within regions similar to those observed to account for its possible radial variations.

### 3.2.2. Model 5Y7M3.7K

For comparison, we report the description of model 5Y7M3.7K, whose initial conditions are significantly different to the ones of 2N7m1K but lead to a similar $f_{enriched}$. In Fig. 7, we show the face-on view of the distribution of the two populations on the left and of the surface density of the whole cluster on the right – both at two evolutionary times – for model 5Y7M3.7K. Comparing it with Fig. 4 of the model 2N7m1K, the central surface
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Density profiles of FG (purple, dashed), SG (orange, dash-dotted), and whole cluster (green, solid) for the model 2N7m1K at three different times, indicated above each panel. The half-mass radii of the three components are reported in the legend. The density profile at 12 Gyr of the whole cluster is reported in black, in the first two panels, for comparison.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Radial profile of mass of SG over the total mass $M_{SG}/M_{tot}$ in black and the fraction of enriched stars within $r$, $f_{enriched}$, in red for the model 2N7m1K at 12Gyr. The vertical grey line represents the half-mass radius of the whole cluster, $r_{h,tot}$.</figcaption></figure>

density is slightly higher here, and, in general, the cluster appears more concentrated at 1 Gyr. At 12 Gyr, the cluster of model 5Y7M3.7K is less concentrated, which is reflected in the larger $r_h$ reported in Table 3. The cluster has lost fewer stars due to both the higher SG initial mass and the larger $W_{0,FG}$. The milder loss of stars can also be seen when looking at the less populated tidal tails.

The profiles of the two models are also quite different; while model 2N7m1K has a very steep profile at the centre after 12 Gyr (see Fig. 5), model 5Y7M3.7K shows a shallower density profile in the inner regions, as shown in Fig. 8. The density profile of the whole cluster is only slightly changed over time, and therefore the final $r_h$ is still very large, not matching the observed values of

Galactic GCs. Also, the total mass of the system is significantly greater than the ones of the bulk of Galactic GCs, meaning that with the initial conditions adopted for this model, the cluster is not losing significant mass but still ends up with an $f_{enriched}$ in the observed range, as shown in Fig. 9. $f_{enriched}$ is well above 0.6 inside 10 pc, while it drops to 0.4 when considering $r > r_{h,tot}$. Here, the difference between the $f_{enriched}$ calculated at the half-mass radius and for the whole cluster is the largest among all models, stressing again the mismatch between the two values.

# 4. Discussion

Through N-body simulations, we explored the effects of various structural and kinematic properties on the mass loss of a massive GC. We now compare our results with the relevant literature and discuss the strengths and limits of our approach.

## 4.1. Comparison with other theoretical works

Understanding whether a cluster can lose a significant number of FG stars, and therefore reproduce the observational constraints discussed in Sect. 1, has been the goal of several studies in the past. Firstly, D’Ercole et al. (2008) addressed the issue, proving the feasibility of such a strong mass loss. It is, however, worth noting that their simulated cluster has an initial mass of $10^4 M_\odot$, significantly smaller than the one we have adopted here, but also of a typical Galactic GC. Clusters with these masses are more prone to lose mass given their shallower potential well; therefore, it is not surprising that they are able to reach higher values of $f_{enriched}$.

Later, Reina-Campos et al. (2018), coupling the Evolution and Assembly of Galaxies and their Environments (EAGLE; Crain et al. 2015) simulations with the subgrid model for stellar cluster formation and evolution MOdelling Star cluster population Assembly In Cosmological Simulations (MOSAICS; Kruijssen et al. 2011), derived that, once assuming an initial $f_{enriched} = 0.05$, their clusters lose a very small number of stars, ending up with fractions of 5–10% for massive systems. Their
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Same as in Fig. 4, but for model 5Y7M3.7K.</figcaption></figure>

result resembles the one we obtained for the model 5Y7m1K, which starts with a similar $f_{\text{enriched}}$ value. In addition, they stated that extrapolating from their results, a present-day massive cluster with $M_{\text{cl}} = 10^6 M_{\odot}$ should initially be composed of a remarkably high fraction of SG stars (~0.8) to reproduce the present-day $f_{\text{enriched}}$ values. This is at variance with what it is generally assumed, and, even in this case, the cluster would be too extended, with a $r_h \sim 10^2$ pc, which is more than one order of magnitude greater than the observed ones. Although in our simulations an extended cluster has to be assumed to match the observed $f_{\text{enriched}}$, the final $r_h$ of the whole cluster reduces significantly during the evolution (in the most extreme case of 2N5m1K, the final cluster $r_h$ is more than 7 times smaller than the FG initial one), which was instead kept fixed in the Reina-Campos et al. (2018) study. On the other hand, we note that a high initial $f_{\text{enriched}}$ does not always lead to a higher final

SG fraction, and therefore we do not need to assume a massive initial SG to reproduce the observed fraction. Such a mismatch between our results and the ones of Reina-Campos et al. (2018) may be ascribed to the positive correlation between the initial and final SG fraction assumed in their work, which we have shown is not always true. In addition, the initial $f_{\text{enriched}}$ may not be the same for clusters of different masses, as was assumed in their study, but a positive trend with cluster mass could be imprinted at birth, as shown by Yaghoobi et al. (2022).

Recently, Sollima et al. (2022) studied the evolution of multiple stellar populations using the binary fraction as a tool to recover the initial concentration required for the SG. Although we do not include the treatment of binaries, model 2N7m1K nicely satisfies the relation required to reproduce the present-day SG binary fraction found by Sollima et al. (2022, their Eq. (3)). Similarly, Vesperini et al. (2021) and Sollima (2021) modelled
<figure><img src="image_9.png" /><figcaption>Fig. 8. Same as in Fig. 5, but for model 5Y7M3.7K.</figcaption></figure>

clusters with masses of the order of $10^6 M_\odot$, reproducing the observed SG fractions, much higher than the ones retrieved in this work (Hypki et al. 2022). For a thorough comparison, other tests should be performed at lower cluster masses.

## 4.2. The fraction of SG lost in the disc

In Table 3, we list the fraction of unbound SG over the total mass of unbound stars, $f_{SG}^{\text{lost}}$. The models can be distinguished into two subgroups: when assuming a massive SG, the unbound stars are composed of ~15% of the SG, while a low-mass SG leads to fractions of ~5%. This quantity is important to study GC evolution in terms of mass loss and its contribution to the Galaxy. Ongoing loss of SG stars was recently identified both in the tails of the disrupting GC Palomar 5 (Phillips et al. 2022) and in the bulge cluster NGC 6723 (Fernández-Trincado et al. 2021). Due to the limited sample, a robust fraction of SG in the tails of these two objects cannot be derived, but it would be useful to add further constraints on the cluster mass loss.

On the other hand, in the Galactic halo field several observational studies have searched for SG-like stars identified by their peculiar chemical composition. In particular, recent investigations have determined fractions of ~2–5% (Carretta et al. 2010; Martell & Grebel 2010; Martell et al. 2011, 2016; Ramírez et al. 2012; Koch et al. 2019). However, this is just a lower limit to the fraction of SG stars that are lost by GCs, due to the high uncertainties regarding the fraction of halo stars that formerly belonged to GCs. Further difficulties arise when focusing on the Galactic disc, where measurements of an SG fraction are not available at present. With the upcoming arrival of 4MOST and WEAVE coupled with Gaia, new insights into the contribution of GCs to the Galactic disc will be provided.

## 4.3. High-z proto-clusters

The mass and half-mass radius of present-day clusters are much smaller than the initial values we assumed for our simulated clusters; however, our assumptions may not be too distant from the real conditions at birth. Accessing the properties of star-forming young GCs is now possible by exploiting gravitational

<figure><img src="image_10.png" /><figcaption>Fig. 9. Same as in Fig. 6, but for model 5Y7M3.7K.</figcaption></figure>

lensing, which makes it possible to identify faint stellar objects at high redshift, in the epoch of their formation (Vanzella et al. 2017). Recently, some proto-GC candidates have been identified, such as the ones in the extended star-forming region strongly magnified by the galaxy cluster MACS J0416.1–2403 (Vanzella et al. 2019; Calura et al. 2021). The region is dominated by two star-forming systems: D1, which has a stellar mass of $2.2 \times 10^7 M_\odot$ and a size of 44 pc, and T1, which is less massive, with its $2 \times 10^6 M_\odot$ and a size of <30 pc. Interestingly, D1 also shows a nucleated star-forming region surrounded by a diffused component. More extended samples of lensed clumps have been presented (Meštrić et al. 2022; Vanzella et al. 2022; Claeyssens et al. 2023) that were detected in various lensed fields and across a wide redshift range (from z ~ 1 to z ~ 8). These samples are composed of clumps of 10–100 pc in size and with masses between $10^5 M_\odot$ and $10^9 M_\odot$, therefore including systems of sizes and masses in the range of our models.
An open problem is determining whether the observed sys-
tems represent single star clusters, extended star-forming com-
plexes, super star clusters (SSCs), or dwarf galaxies. In the MPs
framework, the idea that GCs may form in hierarchical com-
plexes or SSCs is not new (Bekki et al. 2017). Young GCs might
be embedded in a larger structure with similar properties, and
a portion of the parent galaxy or SSCs may provide processed
materials for the creation of MPs (Renzini et al. 2022).

The *James Webb* Space Telescope is opening a new window to the high-redshift observations of GCs. Besides compact clumps and young proto-GCs at high redshift, a recent, exciting discovery has revealed the presence of quiescent, evolved, and massive GCs associated with their host galaxy in the Sparkler system (Mowla et al. 2022). Considering that we are only in the earliest stages of calibration of in-flight data from the JWST, we have exciting times ahead of us, as it is presumable that the current samples may grow rapidly and provide new, fundamental insights into the formation of GCs.

4.4. *Model limitations*

We studied the evolution of a massive cluster orbiting the Milky
Way to explore whether it can lose mass as a result of tidal
effects of the Galactic potential. We derive that, in order to repro-
duce present-day clusters, an initially very extended FG has to
be assumed. However, we have not included ingredients that
are known to increase the mass loss rate, such as gravitational
and tidal shocks, dynamical friction, and the presence of bina-
ries, and dark remnants and do not change the orbital parame-
ters. The addition of these processes would likely increase the
mass loss in our simulations and possibly increase the final
$f_{\text{enriched}}$. On the other hand, the assumed static Galactic potential
overestimates the tidal field acting on the GC, which has been
shown to be much weaker at early times (Renaud et al. 2017).
To overcome these limitations, a study of the dynamical evo-
lution of a GC with MPs is to be performed in a fully cosmo-
logical context. Attempts to study the early formation of GCs
in cosmological simulations are being performed (Kimm et al.
2016; Ma et al. 2020; Li & Gnedin 2019), sometimes with res-
olution high enough to study the feedback of individual stars
(Calura et al. 2022). Although still challenging, it is foreseeable
that in the near future, such tools will also allow us to model
MPs and their long-term dynamical evolution.

5. Conclusions

Most of the scenarios proposed so far for the formation of mul-
tiple stellar populations have to deal with the ‘mass budget’
problem. To overcome it, what it is generally assumed is that
clusters were initially more massive, between five and twenty
times more than they are now (see Bastian & Lardo 2018 and
reference therein). As a consequence, during their evolution they
must lose a significant amount of mass in terms of stars to rec-
oncile with the observational values. We investigated, through a
series of N-body simulations, what the conditions are (if any)
for a massive cluster with an initial mass of $M \sim 10^7 M_\odot$
composed of two different populations to undergo a significant mass
loss during its evolution and end up, after 12 Gyr, with structural
properties in agreement with the present-day GC ones. Our clus-
ter is located in the disc of the Galaxy, and it orbits the centre at
4 kpc. It is therefore evolving under the effect of the tidal field
of the Milky Way. We tested the effects of various parameters
on the mass loss and the fraction of SG stars, $f_{\text{enriched}}$, one of the
strongest observational constraints.

Before performing the simulations with two populations, we
investigated the evolution of single-population clusters. These
results are useful to determine the effects of various param-
eters and for a comparison with the ones obtained with two
populations.

We summarise the main results of the work in the following
points.

- Our best model, 2N7m1K, which starts with $W_{0,FG} = 2$ and a low mass SG of $M_{\text{SG}}^{\text{ini}} = 7 \times 10^5 M_{\odot}$, suffers a strong mass loss, particularly in the FG. It predicts a final total mass of $\sim 4 \times 10^5 M_{\odot}$, in agreement with present-day GCs (Baumgardt & Hilker 2018) with $f_{\text{enriched}}$ at $r_{\text{h,tot}}$ of 0.59, which is slightly lower than the average value for clusters with the same mass but comparable with the ones at the lower edge of the observed interval ($\sim 0.6-0.8$, Milone et al. 2017). On the other hand, the $r_{\text{h,tot}}$ = 8 pc is slightly larger than the ones derived for clusters of similar mass. The FG is reduced by almost 98% of its initial mass and the final mass loss factor is around 20.

- The parameters that affect the mass-loss rate most and, in general, the evolution of the clusters, are the degree of primordial segregation, the FG initial concentration as determined by the initial value of the King dimensionless central potential $W_{0,FG}$, the initial mass of the SG, $M_{\text{SG}}^{\text{ini}}$, and the initial half-mass radius of the SG, $r_{\text{h,SG}}$. In order to lose enough mass, a $W_{0,FG} = 2$ and a low-mass SG of $M_{\text{SG}}^{\text{ini}} = 7 \times 10^5 M_{\odot}$ have to be assumed. Clusters with these initial conditions are able to lose more than 90% of their FG mass, as required to solve the mass budget problem. Such a small $W_{0,FG}$ implies an extended FG with $r_{\text{h}} = 60$ pc, which is comparable to the size of diffuse star clusters observed in high-redshift star-forming complexes (Meštrić et al. 2022; Claeysens et al. 2023).

- From a comparison between the single population models and the two population ones with $W_{0,FG} = 2$, it has been shown that the presence of an SG, even if not very massive, prevents the disruption of the system, as happens when no SG is included.

- Clusters with an initially higher SG mass, and therefore with a higher SG fraction, $f_{\text{enriched}}$, do not always show a higher final SG fraction with respect to clusters starting with a low-mass SG. This is particularly true when small values of $W_{0,FG}$ are assumed. Such behaviour suggests that a positive correlation between the initial and final $f_{\text{enriched}}$ may not always be verified.

- Our clusters are all initially composed of a centrally con-
centrated SG. This difference between the spatial distribu-
tion of the two populations is also found at the end of the
simulations for all models. Consequently, $f_{\text{enriched}}$ is not flat
as a function of radius, and, in particular, it is higher at
the centre and decreases moving outwards. Since observa-
tions are hardly ever able to derive $f_{\text{enriched}}$ for the whole
cluster, but rather for some fraction of $r_h$ only, caution has
to be taken when comparing the simulation results with
observed values. In our cases, differences of up to 20% have
been found between $f_{\text{enriched}}$ of the whole cluster and $f_{\text{enriched}}$
at $r_{\text{h,tot}}$.

- Clusters with a low-mass SG lose a small fraction of SG stars, generally between 4 and 5% of all unbound stars in the tails. On the other hand, clusters with initially massive SG lose 15% of SG stars. These values may be used for comparison with GCs where tidal tails have been detected, such as Palomar 5 and NGC 6723.
Possible follow-ups of the current work could be to expand it, exploring how the intensity of mass loss depends on the initial properties of both the FG and the SG components (e.g., initial masses, galactocentric distance, galaxy potential). A promising method to perform a large series of simulations would be to run a single one-component model and then interpret it a posteri- ori as a multi-component system, a technique recently applied by Nipoti et al. (2021) to a two-component dwarf galaxy orbit- ing the Milky Way (see also Bullock & Johnston 2005; Errani et al. 2015).

Another important step would be to implement a more sophisticated treatment of stellar evolution, exploiting popula- tion synthesis codes such as POSYDON (Andrews 2022) or SEVN (Spera et al. 2019; Iorio et al. 2023) with the latter offer- ing the possibility to track the chemical evolution, which is fun- damental for the study of multiple stellar populations. Moreover, to achieve a more realistic modelling of the phenomenon, the implementation of other physical processes (e.g., binaries, natal kicks, shocks) not taken into account here will be crucial, given that these could potentially trigger more mass loss.

The models could also be adapted to study external galaxies, such as the Magellanic Clouds, where dynamically younger and FG-dominated globular clusters have been found. Such studies will contribute to refining our knowledge of stellar cluster evo- lution and the assembly history of the Galaxy.

Acknowledgements. The authors thank the anonymous referee for a very con- structive report and suggestions that helped significantly improve the quality of the manuscript. EL acknowledges financial support from the European Research Council for the ERC Consolidator grant DEMOBLACK, under contract no. 770017. This work has received funding from INAF Research GTO-Grant Normal RSN2-1.05.12.05.10 – Understanding the formation of globular clusters with their multiple stellar generations (ref. Anna F. Marino) of the “Bando INAF per il Finanziamento della Ricerca Fondamentale 2022”. Part of the calculations presented in this paper were enabled by resources provided by the Swedish National Infrastructure for Computing (SNIC) at Tetrahil and LUNARC. Those resources are partially funded by the Swedish Research Council through grant agreement no. 2018-05973. AMB acknowledges funding from the Euro- pean Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 895174. FC acknowledges support from grant PRIN MIUR 2017-20173ML3WW 001, from the INAF main-stream (1.05.01.86.31) and from PRIN INAF 1.05.01.85.01.

References

Agrawal, P., Hurley, J., Stevenson, S., Szécsi, D., & Flynn, C. 2020, MNRAS, 497, 4549

Allen, C., & Santillan, A. 1991, Rev. Mex. Astron. Astrofis., 22, 255

Ambartsumian, V. A. 1938, TsAGI Uchenye Zapiski, 22, 19

Andrews, J. 2022, AAS/High Energy Astrophysics Division, 54, 110.54

Armandroff, T. E. 1989, AJ, 97, 375

Armandroff, T. E., & Zinn, R. 1988, AJ, 96, 92

Arunima, A., Pfalzner, S., & Govind, A. 2023, A&A, 670, A128

Banerjee, S., & Kroupa, P. 2011, ApJ, 741, L12

Bastian, N., & Lardo, C. 2015, MNRAS, 453, 357

Bastian, N., & Lardo, C. 2018, ARA&A, 56, 83

Bastian, N., Cabrera-Ziri, I., Davies, B., & Larsen, S. S. 2013, MNRAS, 436, 2852

Baumgardt, H., & Hilker, M. 2018, MNRAS, 478, 1520

Baumgardt, H., & Kroupa, P. 2007, MNRAS, 380, 1589

Baumgardt, H., & Makino, J. 2003, MNRAS, 340, 227

Baumgardt, H., De Marchi, G., & Kroupa, P. 2008, ApJ, 685, 247

Bekki, K. 2019, MNRAS, 486, 2570

Bekki, K., Jeřábková, T., & Kroupa, P. 2017, MNRAS, 471, 2242

Bellini, A., Vesperini, E., Piotto, G., et al. 2015, ApJ, 810, L13

Bellini, A., Libralato, M., Bedin, L. R., et al. 2018, ApJ, 853, 86

Bica, E., Bonatto, C., Barbuy, B., & Ortolani, S. 2006, A&A, 450, 105

Bica, E., Ortolani, S., & Barbuy, B. 2016, PASA, 33, e028

Binney, J., & Tremaine, S. 2008, Galactic Dynamics, 2nd edn. (Princeton: Princeton University Press)

Bonetti, M., Bortolas, E., Lupi, A., & Dotti, M. 2021, MNRAS, 502, 3554

Breen, P. G. 2018, MNRAS, 481, L110

Bullock, J. S., & Johnston, K. V. 2005, ApJ, 635, 931

Cabrera-Ziri, I., Bastian, N., Longmore, S. N., et al. 2015, MNRAS, 448, 2224

Calura, F., D'Ercole, A., Vesperini, E., Vanzella, E., & Sollima, A. 2019, MNRAS, 489, 3269

Calura, F., Vanzella, E., Carniani, S., et al. 2021, MNRAS, 500, 3083

Calura, F., Lupi, A., Rosdahl, J., et al. 2022, MNRAS, 516, 5914

Capuzzo-Dolcetta, R., Mastrobuono-Battisti, A., & Maschietti, D. 2011, New Astron., 16, 284

Carretta, E., Bragaglia, A., Gratton, R. G., et al. 2009, A&A, 505, 117

Carretta, E., Bragaglia, A., Gratton, R. G., et al. 2010, A&A, 516, A55

Chernoff, D. F., & Weinberg, M. D. 1990, ApJ, 351, 121

Claeysens, A., Adamo, A., Richard, J., et al. 2023, MNRAS, 520, 2180

Contenta, F., Varri, A. L., & Heggie, D. C. 2015, MNRAS, 449, L100

Cordero, M. J., Pilachowski, C. A., Johnson, C. I., et al. 2014, ApJ, 780, 94

Cordero, M. J., Hénault-Brunet, V., Pilachowski, C. A., et al. 2017, MNRAS, 465, 3515

Cordoni, G., Milone, A. P., Mastrobuono-Battisti, A., et al. 2020, ApJ, 889, 18

Côté, P. 1999, AJ, 118, 406

Crain, R. A., Schaye, J., Bower, R. G., et al. 2015, MNRAS, 450, 1937

Dalessandro, E., Cadelano, M., Vesperini, E., et al. 2019, ApJ, 884, L24

de Mink, S. E., Pols, O. R., Langer, N., & Izzard, R. G. 2009, A&A, 507, L1

Decressin, T., Charbonnel, C., & Meynet, G. 2007, A&A, 475, 859

Denissenkov, P. A., & Hartwick, F. D. A. 2014, MNRAS, 437, L21

D'Ercole, A., Vesperini, E., D'Antona, F., McMillan, S. L. W., & Recchi, S. 2008, MNRAS, 391, 825

D'Ercole, A., D'Antona, F., Ventura, P., Vesperini, E., & McMillan, S. L. W. 2010, MNRAS, 407, 854

D'Ercole, A., D'Antona, F., & Vesperini, E. 2016, MNRAS, 461, 4088

Di Matteo, P. 2016, PASA, 33

Di Matteo, P., Spite, M., Haywood, M., et al. 2020, A&A, 636, A115

Dondoglio, E., Milone, A. P., Lagioia, E. P., et al. 2021, ApJ, 706, 96

D'Orazi, V., Gratton, R., Lucatello, S., et al. 2010, ApJ, 719, L213

Elmegreen, B. G. 2017, ApJ, 836, 80

Errani, R., Penarrubia, J., & Tormen, G. 2015, MNRAS, 449, L46

Fernández-Trincado, J. G., Beers, T. C., Minniti, D., et al. 2021, A&A, 647, A64

Fernández-Trincado, J. G., Beers, T. C., Barbuy, B., et al. 2022, A&A, 663, A126

Fuji, M. S., & Portegies Zwart, S. 2011, Science, 334, 1380

Gieles, M., Charbonnel, C., Krause, M. G. H., et al. 2018, MNRAS, 478, 2461

Giersz, M., Askar, A., Wang, L., et al. 2019, MNRAS, 487, 2412

Gnedin, O. Y., & Ostriker, J. P. 1997, ApJ, 474, 223

Gratton, R., Bragaglia, A., Carretta, E., et al. 2019, A&ARv, 27, 8

Haghi, H., Hoseini-Rad, S. M., Zonoozi, A. H., & Küpper, A. H. W. 2014, MNRAS, 444, 3699

Harris, W. E. 2010, arXiv e-prints [arXiv:1012.3224]

Hénault-Brunet, V., Gieles, M., Agertz, O., & Read, J. I. 2015, MNRAS, 450, 1164

Hypki, A., Giersz, M., Hong, J., et al. 2022, MNRAS, 517, 4768

Iorio, G., Mapelli, M., Costa, G., et al. 2023, MNRAS, 524, 426

Jang, S., Milone, A. P., Legnardi, M. V., et al. 2022, MNRAS, 517, 5687

Kamann, S., Dalessandro, E., Bastian, N., et al. 2020, MNRAS, 492, 966

Khalaj, P., & Baumgardt, H. 2015, MNRAS, 452, 924

Khalaj, P., & Baumgardt, H. 2016, MNRAS, 457, 479

Kimm, T., Cen, R., Rosdahl, J., & Yi, S. K. 2016, ApJ, 823, 52

King, I. R. 1966, AJ, 71, 64

Koch, A., Grebel, E. K., & Martell, S. L. 2019, A&A, 625, A75

Kroupa, P. 2001, MNRAS, 322, 231

Kruijssen, J. M. D. 2015, MNRAS, 454, 1658

Kruijssen, J. M. D., Pelupessy, F. I., Lamers, H. J. G. L. M., Portegies Zwart, S. F., & Icke, V. 2011, MNRAS, 414, 1339

Krumholz, M. R., McKee, C. F., & Bland-Hawthorn, J. 2019, ARA&A, 57, 227

Küpper, A. H. W., Kroupa, P., & Baumgardt, H. 2008, MNRAS, 389, 889

Küpper, A. H. W., Maschberger, T., Kroupa, P., & Baumgardt, H. 2011, MNRAS, 417, 2300

Lacchin, E., Calura, F., & Vesperini, E. 2021, MNRAS, 506, 5951

Lacchin, E., Calura, F., Vesperini, E., & Mastrobuono-Battisti, A. 2022, MNRAS, 517, 1171

Lada, C. J., Margulis, M., & Dearborn, D. 1984, ApJ, 285, 141

Lardo, C., Bellazzini, M., Pancino, E., et al. 2011, A&A, 525, A114

Larsen, S. S., Strader, J., & Brodie, J. P. 2012, A&A, 544, L14

Larsen, S. S., Brodie, J. P., Grundahl, F., & Strader, J. 2014, ApJ, 797, 15

Lee, J.-W. 2015, ApJS, 219, 7

Lee, J.-W. 2017, ApJ, 844, 77

Lee, J.-W. 2018, ApJS, 238, 24

Leigh, N. W. C., Mastrobuono-Battisti, A., Perets, H. B., & Böker, T. 2014, MNRAS, 441, 919

Li, H., & Gnedin, O. Y. 2019, MNRAS, 486, 4030

Libralato, M., Bellini, A., Piotto, G., et al. 2019, ApJ, 873, 109

Libralato, M., Vesperini, E., Bellini, A., et al. 2023, ApJ, 944, 58
Lucatello, S., Sollima, A., Gratton, R., et al. 2015, A&A, 584, A52

Ma, X., Grudić, M. Y., Quataert, E., et al. 2020, MNRAS, 493, 4315

Marino, A. F., Milone, A. P., Renzini, A., et al. 2019, MNRAS, 487, 3815

Martell, S. L., & Grebel, E. K. 2010, A&A, 519, A14

Martell, S. L., Smolinski, J. P., Beers, T. C., & Grebel, E. K. 2011, A&A, 534, A136

Martell, S. L., Shetrone, M. D., Lucatello, S., et al. 2016, ApJ, 825, 146

Masseron, T., García-Hernández, D. A., Mészáros, S., et al. 2019, A&A, 622, A191

Mastrobuono-Battisti, A., Di Matteo, P., Montuori, M., & Haywood, M. 2012, A&A, 546, L7

Mastrobuono-Battisti, A., Khoperskov, S., Di Matteo, P., & Haywood, M. 2019, A&A, 622, A86

McLaughlin, D. E., & van der Marel, R. P. 2005, ApJS, 161, 304

Meštrić, U., Vanzella, E., Zanella, A., et al. 2022, MNRAS, 516, 3532

Milone, A. P., & Marino, A. F. 2022, Universe, 8, 359

Milone, A. P., Piotto, G., Bedin, L. R., et al. 2012, ApJ, 744, 58

Milone, A. P., Piotto, G., Renzini, A., et al. 2017, MNRAS, 464, 3636

Milone, A. P., Vesperini, E., Marino, A. F., et al. 2020, MNRAS, 492, 5457

Minniti, D. 1995, AJ, 109, 1663

Mowla, L., Iyer, K. G., Desprez, G., et al. 2022, ApJ, 937, L35

Nipoti, C., Cherchi, G., Iorio, G., & Calura, F. 2021, MNRAS, 503, 4221

Norris, J., & Freeman, K. C. 1979, ApJ, 230, L179

Phillips, S. G., Schiavon, R. P., Mackereth, J. T., et al. 2022, MNRAS, 510, 3727

Piotto, G., Villanova, S., Bedin, L. R., et al. 2005, ApJ, 621, 777

Pouliasis, E., Di Matteo, P., & Haywood, M. 2017, A&A, 598, A66

Ramírez, I., Meléndez, J., & Chanamé, J. 2012, ApJ, 757, 164

Reina-Campos, M., Kruijssen, J. M. D., Pfeffer, J., Bastian, N., & Crain, R. A. 2018, MNRAS, 481, 2851

Renaud, F., Agertz, O., & Gieles, M. 2017, MNRAS, 465, 3622

Renzini, A., D'Antona, F., Cassisi, S., et al. 2015, MNRAS, 454, 4197

Renzini, A., Marino, A. F., & Milone, A. P. 2022, MNRAS, 513, 2111

Richer, H. B., Heyl, J., Anderson, J., et al. 2013, ApJ, 771, L15

Schaerer, D., & Charbonnel, C. 2011, MNRAS, 413, 2297

Schiavon, R. P., Zamora, O., Carrera, R., et al. 2017, MNRAS, 465, 501

Simioni, M., Milone, A. P., Bedin, L. R., et al. 2016, MNRAS, 463, 449

Sollima, A. 2021, MNRAS, 502, 1974

Sollima, A., Ferraro, F. R., Bellazzini, M., et al. 2007, ApJ, 654, 915

Sollima, A., Nipoti, C., Mastrobuono Battisti, A., Montuori, M., & Capuzzo-Dolcetta, R. 2012, ApJ, 744, 196

Sollima, A., Gratton, R., Lucatello, S., & Carretta, E. 2022, MNRAS, 512, 776

Spera, M., Mapelli, M., Giacobbo, N., et al. 2019, MNRAS, 485, 889

Spitzer, L., & J., 1940, MNRAS, 100, 396

Szigeti, L., Mészáros, S., Szabó, G. M., et al. 2021, MNRAS, 504, 1144

Tanikawa, A., & Fukushige, T. 2009, PASJ, 61, 721

Teuben, P. 1995, in Astronomical Data Analysis Software and Systems IV, eds. R. A. Shaw, H. E. Payne, & J. J. E. Hayes, ASP Conf. Ser., 77, 398

Van Den Bergh, S. 2003, ApJ, 590, 797

Vanzella, E., Calura, F., Meneghetti, M., et al. 2017, MNRAS, 467, 4304

Vanzella, E., Calura, F., Meneghetti, M., et al. 2019, MNRAS, 483, 3618

Vanzella, E., Castellano, M., Bergamini, P., et al. 2022, ApJ, 940, L53

Vesperini, E., & Heggie, D. C. 1997, MNRAS, 289, 898

Vesperini, E., McMillan, S. L. W., & Portegies Zwart, S. 2009, ApJ, 698, 615

Vesperini, E., McMillan, S. L. W., D'Antona, F., & D'Ercole, A. 2010, ApJ, 718, L112

Vesperini, E., McMillan, S. L. W., D'Antona, F., & D'Ercole, A. 2013, MNRAS, 429, 1913

Vesperini, E., Hong, J., Giersz, M., & Hypki, A. 2021, MNRAS, 502, 4290

Vincenzo, F., Matteucci, F., Belfiore, F., & Maiolino, R. 2016, MNRAS, 455, 4183

von Hoerner, S. 1957, ApJ, 125, 451

Wang, L., Kroupa, P., Takahashi, K., & Jerabkova, T. 2020, MNRAS, 491, 440

Webb, J. J., Harris, W. E., Sills, A., & Hurley, J. R. 2013, ApJ, 764, 124

Yaghoobi, A., Calura, F., Rosdahl, J., & Haghi, H. 2022, MNRAS, 510, 4330

Zennaro, M., Milone, A. P., Marino, A. F., et al. 2019, MNRAS, 487, 3239

Zinn, R. 1985, ApJ, 293, 424