# Simulations of Globular Clusters Within Their Parent Galaxies: Multiple
Stellar Populations And Internal Kinematics

Madeleine McKenzie,1 Kenji Bekki,1
1 ICRAR, M468, The University of Western Australia, 35 Stirling Highway,
Crawley, WA 6009, Australia E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Using three-dimensional smoothed particle hydrodynamics simulations, we
investigate the formation of multiple stellar populations (MSPs) in globular
clusters (GCs) within the context of their parent galaxies. In our scenario,
the second generation (2G) of stars originate from both asymptotic giant
branch (AGB) polluters and pristine gas accreted from the host galaxy.
Previous theoretical and numerical studies have demonstrated that this “AGB
with dilution” model has the potential to alleviate several problems faced by
the classical AGB scenario. However, the accretion of pristine gas onto the GC
has yet to be investigated within the context of the parent galaxy. This paper
presents the preliminary results from our original simulation code which
models GC formation from giant molecular clouds in a host galaxy, and
subsequent gas accretion onto the GC. By simulating the genesis of the 2G over
a 370 Myr time frame, we demonstrate that the fraction of 2G stars is
inextricably linked to the GC’s environment. Our simulations rationalise the
wide variety of abundance patterns, kinematics and 2G concentrations by
altering the initial conditions of both the GC progenitor and the host galaxy
itself. Most notably, we reproduce a positive correlation between the fraction
of 2G stars and the initial mass of the cluster. We discuss the physical
implications of our scenario and compare our simulations with observations of
the Galactic GC 47 Tucanae (47 Tuc). Finally, we present scaling relations
which encompass the wider GC population and serve as a reference for future
observations.

###### keywords:

globular clusters: general – hydrodynamics: methods – stars: formation –
galaxies: star clusters: general.

††pubyear: 2020††pagerange: Simulations of Globular Clusters Within Their
Parent Galaxies: Multiple Stellar Populations And Internal Kinematics–B

##  1 Introduction

The answer to how multiple stellar populations (MSPs) form within globular
clusters has evaded astronomers for several decades. Since the suggestion that
MSP may be responsible for the Cyanogen bimodality in 47 Tucanae (47 Tuc) by
Norris & Freeman (1979), a wealth of photometric (e.g. Sbordone et al. 2011;
Piotto et al. 2015; Milone et al. 2017; Lagioia et al. 2019; Lee & Sneden
2020) and spectroscopic (e.g. Barbuy et al. 2009; Carretta et al. 2009a;
Kamann et al. 2018; Latour et al. 2018; Husser et al. 2020) evidence indicates
that MSPs are the source of light element inhomogenities almost ubiquitously
observed in GCs. Several groups and individuals have postulated various
scenarios which fulfil observational criteria to varying degrees, but none
have been universally accepted.

The majority of observed GCs exhibit anti-correlations between Na and O (e.g.
Carretta 2019) and occasionally Al and Mg (Pancino et al. 2017). Subtle
differences in light element abundances are used to delineate different
populations of stars. Enriched populations are characterised by enhanced N, Na
and Al and depleted C, O and Mg. The terminology for such population varies
widely in the literature, with first and second generations (1G and 2G
respectively, introduced in Piotto et al. (2015)) commonly being employed.
Such nomenclature assumes a linear formation scenario which has yet to be
confirmed, however, we adopt it for continuity. 1G stars show primordial
(oxygen-rich, sodium-poor) chemical composition similar to that of field stars
(e.g. Martell 2011). Elemental abundance patterns of 2G stars are ascribed by
CNO cycling and p-capture processes at high temperatures (e.g. Gratton et al.
2012). It has been hypothesised that the chemical abundance patterns of 2G
stars are specific to GCs, however, N-enriched stars akin to the 2G have been
identified by the APOGEE survey within the galactic bulge (Schiavon et al.,
2016). One explanation for this is that these 2G stars are the remains of
tidal interactions between GCs and the Milky Way (Vesperini et al. 2010;
Martell et al. 2016).

Further classifications exist for the general GC population. Type I GCs
constitute $\approx 83\%$ of all clusters (Marino et al., 2019) and show a
narrow spread in iron (e.g. Carretta et al. 2009b; Bailin 2019). Type II GCs
show evidence of several groups of stars with varying iron abundances with
notable examples including $\omega$ Centauri (Lee et al. 1999; Bedin et al.
2004), M54 (Carretta et al. 2010b; Layden & Sarajedini 2000) and NGC1851
(Carretta et al., 2011). Popularised by Milone et al. (2017), pseudo two-color
diagrams known as ‘chromosome maps’ have become a powerful tool for diagnosing
and quantifying multiple populations, for example, four populations have been
identified in NGC2080 (Milone et al. 2015; Latour et al. 2019) and five in M15
(Nardiello et al., 2018) using this technique. By exploiting small chemical
differences between stars, it is clear that the notion of single stellar
populations must be discarded.

Due to the wide variety of physical and dynamical conditions observed in GCs,
a robust yet flexible scenario is required in order to encapsulate their
formation processes. Milone (2019) concisely summarises the observational
criteria for which a scenario must explain. For example, the scenario must be
capable of reproducing a central concentration of 2G stars (e.g. 47 Tuc in
Milone et al. 2018 and NGC3201 in Carretta et al. 2010a), equivalent
concentrations of 1G and 2G stars (e.g. M71 in Gerber et al. 2020) and on
rarer occasions, a centrally concentrated 1G (i.e. M15 in Larsen et al. 2015
or M80 in Dalessandro et al. 2018). It should reproduce a fraction of 2G stars
which varies from as low as 25% (NGC339 in Niederhofer et al. (2017)) to 90%
in $\omega$ Centauri (Bellini et al., 2009). Additionally, there should be a
way of explaining the stellar halo around M2 (Kuzma et al., 2016) and NGC1851
(Kuzma et al., 2018).

Several reviews by Gratton et al. (2012), Renzini et al. (2015), Bastian &
Lardo (2018) and most recently Gratton et al. (2019) summarise the leading
formation scenarios for GCs. Many of the canonical formation scenarios have
been plagued by a variety of unexplained observational patterns. For example,
the AGB scenario first discussed by Iben (1975), Iben (1976) and later in
Cottrell & Da Costa (1981) succeeds in justifying a number of observations
(e.g. Ventura & D’Antona 2008b; Ventura & D’Antona 2008a; Choi & Yi 2008;
D’Orazi & Marino 2010; D’Ercole et al. 2012; D’Antona et al. 2016) but suffers
from two main criticisms. The "mass budget problem" (e.g. Smith & Norris 1982,
Renzini 2008) is where the progenitor mass of the 1G is not large enough to
form a massive 2G. This issue could be alleviated by employing a top-heavy IMF
(Prantzos & Charbonnel 2006; Bekki et al. 2017; Kroupa 2019) or through
assuming a mechanism by which a significant number of primordial stars are
stripped (e.g. Decressin et al. 2008). Secondly, AGB models predict a Na-O
correlation rather than anti-correlation (Bastian & Lardo, 2018). A solution
to these problems, suggested by D’Ercole et al. (2008), was to allow AGB
ejecta to mix with the pristine gas with comparable metallicities to the 1G.
Additional mass is obtained through accretion of gas from the parent galaxy
and the enrichment patterns of 2G stars are thus dependent on the degree of
mass accretion as well as AGB yields. Several other polluters besides AGB
stars have been accused of producing MSPs; including fast-rotating massive
main sequence stars (Decressin et al. 2007), binary mergers (de Mink et al.
2009), supermassive main sequence stars (Denissenkov & Hartwick 2014) and cool
supergiant stars (Szécsi & Wünsch 2019).

Recently, the scenario proposed by Wang et al. (2020) assumes a multichannel
scenario where not one particular polluter is responsible for the formation of
enriched populations. Through dynamical modelling of binary mergers,
supermassive and fast-rotating massive main sequence stars, they generate a
chemically enriched population analogous to observed Galactic GCs. As there is
no need for gas accretion, they state that their scenario is independent of
the galactic conditions of their parent galaxy. However, their scenario is at
odds with recent results from Milone et al. (2020) where they found that
Magellanic Cloud GCs had a notably lower 2G fraction than what was predicted
from the mass to enriched fraction found in Galactic GCs.

The hierarchical formation model of $\Lambda$CDM dictates that our Galaxy
underwent several merging events in its history (e.g. Belokurov et al. 2018;
Helmi et al. 2018; Myeong et al. 2019; Kruijssen et al. 2020). Merging
episodes brought a considerable number of GCs into our galaxy (e.g. Myeong et
al. 2018; Horta et al. 2020; Forbes 2020) several of which have since been
traced back to these events. These GCs have been identified through velocity
and action space calculations, and APOGEE surveys have chemically linked these
clusters as well. If the fraction of enriched stars formed within GCs is
mediated by the host galaxy as found by Milone et al. (2020), a formation
scenario should be capable of accounting for these differences.

Recent hydrodynamical simulations of GC formation have demonstrated that
feedback effects from SNe and AGB stars within fractal giant molecular clouds
(GMCs) influence the chemical and dynamical properties of 2G stars (Bekki
2017a; Bekki 2017b). Calura et al. (2019) recently investigated the mixing
process of accreting interstellar medium (ISM) and AGB ejecta in a self-
consistent manner and reproduced the observed He abundance spread in GCs.
Using new hydrodynamical simulations with a mass resolution of $\sim 0.01{\rm
M}_{\odot}$, Bekki (2019) investigated the impact of direct interaction of
individual 1G stars with intra-cluster gas. Due to computational complexities,
simulations of MSP formation tend to neglect the impact of the parent galaxies
dynamics. This oversimplification hinders the investigation of a number of key
observations such as the correlation between GC masses and 2G fraction (Milone
et al. 2017 and references therein) and internal rotation of the 1G and 2G
populations (e.g. Bianchini et al. 2013 Bellini et al. 2015; Bianchini et al.
2018; Bellini et al. 2018, Cordoni et al. 2020).

The purpose of this paper is to discuss our preliminary findings from our
original simulation code. We present a simulated analogue of the Galactic GC
47 Tuc, an archetypal Type I GC, to illustrate the success and drawbacks of
this model. We trace the evolution of our fiducial 47 Tuc analogue over a 370
Myr window and analyse the morphological, chemical and kinematic consequences
of gas accretion onto the cluster. Next, we extend our analysis to a broader
population of GC and parent galaxy combinations and demonstrate that the gas
fraction and density of the galaxy mediates the relationship between the
fraction of 2G stars and total mass of the cluster as found by Milone et al.
(2020). Given the presence of rotation in a number of clusters (e.g. see
Sollima et al. 2019) we describe the correlations between the 1G and 2G
kinematics and the total cluster mass. Finally, we discuss how these
simulations can naturally give rise to multiple discrete 2G populations and
Type II GCs.

The simulation code employed in the present study allows us to investigate how
complex stellar clusters can form within their natal GMCs in dwarf galaxies.
Previous theoretical and numerical studies have investigated the physical
roles of external gas accretion onto GCs in detail (e.g., Bekki & Mackey 2009;
Pflamm-Altenburg & Kroupa 2009; Conroy & Spergel 2011; Armstrong et al. 2018;
McKenzie & Bekki 2018). However, these works did not self-consistently
investigate the formation of 1G stars and subsequent gas accretion on the
existing stellar systems. Thus it is unclear how additional populations of
stars can form and how the properties of the parent galaxies impact such
populations. Over a series of papers, we intend to develop our understanding
of GC formation and trace the impact of the parent galaxy throughout the
cluster’s evolution.

We structure the paper as follows. In Section 2, we describe our original
simulation code. In Section 3, we present a 47 Tuc analogue for our scenario
and global scaling relations from our scenario. In Section 4, we discuss the
ramifications of our simulations and summarise our conclusions in Section 5.

##  2 The model

###  2.1 A GC formation scenario

In the standard AGB scenario, an essential physical process of GC formation is
the efficient conversion of ISM mixed with AGB ejecta into new stars. In the
present study, we assume these mixing processes and subsequent star formation
can occur within high redshift gas-rich dwarf galaxies. We use galaxy-scale
hydrodynamical simulations to investigate the formation of 1G and 2G stars
simultaneously; however, some physical processes are poorly understood (e.g.,
dynamical and thermal influences of multiple SNe on GC-hosting GMCs).
Accordingly, we adopt a somewhat idealized model for various physical
processes related to GC formation. We describe these simplifications and
discuss possible shortcomings of our GC formation model in Section 4 and
suggest ways of constructing more realistic models in future works.

####  2.1.1 Formation of massive molecular clouds in dwarfs

In order to explain the observed typical mass ($\sim 3\times 10^{5}{\rm
M}_{\odot}$) and high fraction of 2G stars ($\sim 0.7$), the original mass of
GCs has been suggested to be $1-3\times 10^{6}{\rm M}_{\odot}$ (the so-called
mass budget problem). This implies that natal GMCs should be at least
$[5-15]\times 10^{6}{\rm M}_{\odot}$ when assuming a star formation efficiency
within the GMCs of around 20%. Therefore, only gas-rich dwarf galaxies with
the potential to form massive GMCs can host GCs in the present formation
scenario. We attribute the formation of massive GMCs to gravitational
instabilities in the stellar and gaseous discs of dwarf galaxies with high
baryonic fractions $f_{\rm b}$. Hierarchical merging of numerous smaller GMCs
is required to generate a cloud more massive than $10^{6}{\rm M}_{\odot}$
which may eventuate into a GC. Merging or collisions of two GMCs may trigger
the formation of massive OB stars in nearby galaxies (e.g., Fukui et al. 2017;
Tsuge et al. 2019). Hence, multiple merging or collision events could produce
a large number of OB stars which may influence the evolution of GC-hosting
GMCs through energetic feedback effects.

####  2.1.2 Formation of 1G stellar systems and the subsequent dispersion of
GMCs by 1G SNe

Following the creation of massive GMCs, 1G stars within a cloud can
efficiently form within short time-scales of a few Myrs. However, the
formation process leaves behind a significant fraction of the original gas. 1G
stars quickly develop into compact stellar systems, although the dynamical
relaxation process is still incomplete after 10 Myr. The most massive of these
stars ($m\sim[50-100]{\rm M}_{\odot}$) explode as Type II SNe in $\sim 3$ Myr
after the formation of 1G stars. A large amount of thermal and kinetic energy
from these multiple SNe can disperse the remaining gas completely, and the SN
feedback effects can form a giant gaseous hole in its host dwarf galaxy. Such
giant HI holes are observed in some dwarf galaxies, and this is discussed
further in Section 4.4.

####  2.1.3 Accretion of ISM and AGB ejecta onto 1G systems

After the complete dispersal of GC-hosting GMCs, 1G stellar systems can
gravitationally attract their surrounding ISM and provided that they are
massive enough, can accrete a significant portion of gas. The stellar winds of
massive AGB stars can also be gravitationally trapped by the stellar systems
and can mix with the accreted ISM. The accretion process of ISM can proceed
sporadically, and it is not like “Bondi-type” (Bondi, 1952) continuous
accretion as the GC-host dwarfs contain clumpy ISM consisting of numerous
small and large gas clouds. This clumpy accretion can result in discrete
multiple stellar populations, as discussed later.

####  2.1.4 2G formation within 1G systems

Once the 1G accumulates some critical mass of ISM and AGB ejecta in its
central region, new stars can be formed efficiently from the gas within the
gravitational potentials of the systems. These 2G stars can have more compact
spatial distributions compared to their 1G counterparts, as demonstrated in
previous numerical simulations of GC formation with MSPs (e.g., Bekki 2011).
However, it could be possible that 2G stars can have less compact
distributions than the 1G provided that the accreting ISM has a considerable
amount of angular momentum with respect to the centre of 1G stellar systems.
The mass-ratios of 2G to 1G stars can possibly depend on the ISM accretion
processes and thus on the physical properties of GC-hosting dwarfs.

###  2.2 Simulation code

In order to investigate the details of the above key formation processes of
GCs with MSPs, we use our original chemodynamical simulation code based on the
SPH method that can run on GPU machines (Bekki 2013; Bekki 2015). This code
allows us to investigate the formation of dust from SNe and AGB stars, the
destruction of dust by SNe in ISM, the formation of molecular hydrogen on dust
grains, chemical evolution of various elements, star formation (on galaxy-
scale), and various feedback effects from stars. However, simulating these
complicated processes of ISM is numerically costly. Because the main focus of
this study is not the formation and evolution of dust and molecular hydrogen,
we switch off calculations related to dust and molecular hydrogen formation.
We also newly introduce our original multiple gravitational lengths with very
narrow maximum time step width to investigate gas accretion processes on GCs
as described below.

###  2.3 Gas-rich dwarf galaxy

A GC-forming dwarf galaxy is assumed to consist of a dark matter halo (DM),
stellar disc, and gaseous disc. The initial total masses of DM, stellar disc,
and gas disc are denoted as $M_{\rm dm}$, $M_{\rm s}$, and $M_{\rm g}$,
respectively. We mainly investigate low-mass dwarf disc galaxies in which (i)
$M_{\rm dm}=[1-5]\times 10^{10}{\rm M}_{\odot}$ and (ii) the mass ratio of the
disc ($M_{\rm s}+M_{\rm g}$) to the dark matter halo ($M_{\rm h}$) in a dwarf
disc galaxy ranges from 0.01 to 0.06. The adopted small baryonic fractions are
quite reasonable, given that recent observations have found that low-mass
halos have smaller $f_{\rm b}$.

We adopt the ‘NFW’ profile for the dark matter halo (Navarro et al. 1996) with
a central cusp predicted by the Cold Dark Matter (CDM) model:

${\rho}(r)=\frac{\rho_{0}}{(r/r_{\rm s})(1+r/r_{\rm s})^{2}},$ (1)

where $r$, $\rho_{0}$, and $r_{\rm s}$ are the distance from the centre of the
cluster, the central density, and the scale-length of the dark halo,
respectively. The virial radius ($r_{\rm vir}$), the scale radius ($r_{\rm
s}$), and the ‘$c$’ parameter (=$r_{\rm vir}/r_{\rm s}$) are chosen such that
the values are consistent with recent cosmological simulations for the adopted
$M_{\rm h}$ (Neto et al. 2007). We mainly investigate the models with $c=16$,
which is reasonable for low-mass halos.

The mass and size of the galactic bulge in a disc galaxy are free parameters
denoted as $M_{\rm b}$ and $R_{\rm b}$, respectively. The radial ($R$) and
vertical ($Z$) density profiles of the adopted exponential stellar disc are
assumed to be proportional to $\exp(-R/R_{0})$ with scale length
$R_{0}=0.2R_{\rm s}$ and to ${\rm sech}^{2}(Z/Z_{0})$ with scale length
$Z_{0}=0.04R_{\rm s}$, respectively. The gas disc with a size $R_{\rm
g}=R_{\rm s}$ has the radial and vertical scale lengths of $0.2R_{\rm g}$ and
$0.02R_{\rm g}$, respectively. The disc of the present model has $R_{\rm
s}=17.5$ kpc and in addition to the rotational velocity caused by the
gravitational field of disc, bulge, and dark halo components. The initial
radial and azimuthal velocity dispersions are assigned to the disc component
according to the epicyclic theory with Toomre’s parameter $Q$ = 1.5. The gas
mass fraction ($f_{\rm g}=M_{\rm g}/M_{\rm s}$) is also a free parameter in
the present study.

We mainly investigate models where the initial gaseous temperature and
metallicity in the ISM is set to 100K and [Fe/H]=$-0.8$ respectively, as we
wish to find the most suitable model for 47 Tuc with the observed metallicity
of [Fe/H]$\sim-0.8$. However, we investigate the models with different [Fe/H]
in the ISM. The radiative cooling processes are properly included by using the
cooling curve by Rosen & Bregman (1995) for $T<10^{4}$K and the MAPPING III
code for $T\geq 10^{4}$K (Sutherland & Dopita, 1993). Heating of the ISM by
dust may be significant in its evolution within galaxies (e.g., photo-electric
heating; Osman et al. 2020), but has been excluded due to the cost incurred by
additional calculations in the hydrodynamical simulations.

###  2.4 Structure and kinematics of the 1G system

The 1G stellar system of the GC is assumed to have a mass of $M_{\rm gc}$, a
size of $R_{\rm gc}$, and a radial scale-length of $0.2R_{\rm gc}$. The GC is
represented by a rotating Plummer model where $M_{\rm gc}$ and $R_{\rm gc}$
are free parameters. The GC is assumed to be supported by velocity dispersion
and rotation and the dispersion is assumed to be isotropic. We introduce the
following parameter to describes the ratio of total rotational energy ($T_{\rm
rot}$) of a GC to its total kinetic energy ($T_{\rm kin}$):

$f_{\rm rot}=\frac{T_{\rm rot}}{T_{\rm kin}}.$ (2)

We assume that the GC’s spin axis is inclined by $\theta_{\rm gc}$ degrees
with respect to the spin axis of its host galaxy (the $z$-axis).

###  2.5 Giant HI hole after GC formation

A very high star formation efficiency within the GMCs is required to allow for
the creation of bound star clusters. A small fraction of GC-forming gas clouds
may remain within the cluster after 1G formation; however, multiple SNe from
the 1G should expel all of the remaining gas. Thus we assume that the natal
GMC of the cluster is both consumed by efficient 1G formation and dispersed by
energetic SN feedback effects after 1G formation. Our current hydrodynamical
simulation can not achieve the required spatial resolution to investigate how
multiple SNe events starting from just 3 Myr after the 1Gs creation can
influence the remaining gas trapped within the cluster. We therefore
approximate this process by artificially inducing a giant gaseous hole after
1G formation centred on the birthplace of a GC. We admit that this is an over-
simplification of GC formation, however, we believe that this allows us to
investigate how ISM and AGB ejecta can be accreted onto existing 1G systems.
The size of the gas hole ($R_{\rm HI}$) is a free parameter ranging from 100
pc to 1 kpc.

###  2.6 Formation model for 2G stars

In our scenario, new stars can form through the accretion of ISM onto the 1G
or from gas ejected from AGB stars, provided that the gas is gravitationally
trapped within the cluster’s potential. We assume that gas particles can be
converted into collisionless new stellar particles (‘new stars’) if the
following three physical conditions are met. Firstly, the local dynamical
timescale should be shorter than the sound crossing time for a gas particle:
this can mimic the Jeans instability of gas. Secondly, the local density
($\rho_{\rm g}$) around a gas particle must exceed a threshold density
(${\rho}_{\rm th}$) for star formation:

${\rho}_{\rm g}>{\rho}_{\rm th}.$ (3)

We assume that this threshold gas density $\rho_{\rm th}$ is $10^{3}$ H atoms
cm-3 for all models in the present study. Finally, the local velocity field
around a gas particle is consistent with that for gravitationally collapsing,
which we model as follows:

$div{\bf v}<0.$ (4)

###  2.7 Gas ejection from AGB stars

We adopt the same model used in Bekki (2017b) in order to investigate how AGB
ejecta can be accreted onto 1G stellar system. In the adopted model, each AGB
star ejects gas particles with chemical abundances predicted from recent AGB
models (e.g., Karakas (2010)). The ejection of new particles from AGB stars
(‘AGB particle’) in a simulation means that the total number of gas particles
(AGB and ISM particles) can slowly increase as a simulation progresses (after
the 1G formation). We consider the ejection of AGB particles at different
epochs ($t_{\rm agb}$ depending on the progenitor masses of AGB stars. For
example, the five epochs of 200, 120, 80, 60, and 40 Myr correspond to the
lifetimes of stars with masses, (i) $3\leq m_{\rm s}<4$ (${\rm M}_{\odot}$),
(ii) $4\leq m_{\rm s}<5$ (${\rm M}_{\odot}$), (iii) $5\leq m_{\rm s}<6$ (${\rm
M}_{\odot}$), (iv) $6\leq m_{\rm s}<7$ (${\rm M}_{\odot}$), and (v) $7\leq
m_{\rm s}<8$ (${\rm M}_{\odot}$), respectively. We set the minimum mass of AGB
stars to $3{\rm M}_{\odot}$ for three reasons; elemental abundences observed
in Galactic GCs require a small age gap between 1G and 2G stars, the fraction
of gaseous ejecta from AGB stars with $m_{\rm s}<3{\rm M}_{\odot}$ is
neglidgible and we simulate a short time frame of $<1$ Gyr of evolution.

At the end of a stellar particle’s main sequence phase, it ejects an AGB
particle with a wind velocity of $v_{\rm wind}=10$ km s-1. The adopted $v_{\rm
wind}=10$ km s-1 is reasonably consistent with recent observations of AGB
stars in the LMC (e.g., Marshall et al. 2004). The initial temperature of AGB
wind ($T_{\rm wind}$) is set to be 1000 K, which is also consistent with
standard theoretical models of AGB winds. In low-mass GCs, the initially warm
AGB ejecta cannot be converted into new stars owing to the shallow
gravitational potentials.

We do not intend on examining the chemical abundances of 1G and 2G stars in
the present study for three reasons. Our current simulations do not
comprehensively model the ISM abundances within the parent galaxy, and thus in
our next paper, we intend on performing an extensive, self-consistent
investigation into the initial gas abundances of the GC’s host. Secondly,
there are some uncertainties in the nucleosynthesis yields of AGB stars in
different groups (e.g., Karakas & Lattanzio 2014), which may significantly
influence the CNO abundance patterns of 2G stars. Finally, the present study
already contains several new results of GC formation, which we wish to focus
exclusively on in this paper.

###  2.8 Multiple gravitational softening lengths

In order to investigate gas accretion processes onto GCs within a galactic-
scale simulation, we adopt vastly different gravitational softening lengths
for dark matter ($\epsilon_{\rm dm}$), stellar and gaseous discs
($\epsilon_{\rm s}$), GCs ($\epsilon_{\rm gc}$). It is challenging to
investigate the long-term dynamical evolution driven by two-body relaxation
for massive GCs with $M_{\rm gc}\sim 10^{6}{\rm M}_{\odot}$ (1G systems
consisting of more than $10^{6}$ stars) in live galactic gravitational
potentials using NBody6 codes, even if there is no gas dynamics (e.g., Rossi
et al. 2016). Accordingly, we introduce a rather small $\epsilon_{\rm gc}$ for
GCs to investigate only the short-term ($\sim 10^{8}$ yr) gas accretion
processes. The long-term evolution of 1G and 2G stars has been demonstrated to
be a key process that governs the spatial structures and internal kinematics
of GCs with MSPs (e.g., Vesperini et al. 2018). Clearly, these important long-
term dynamical evolution of 1G and 2G stars in GCs should be investigated in
future studies with a proper simulation code (e.g., Nbody6).

###  2.9 Parameter study

We mainly describe the results of our fiducial model (M1) with $M_{\rm
dm}=5\times 10^{9}\ {\rm M}_{\odot}$, $R_{\rm vir}=24$ kpc, $f_{\rm b}=0.036$,
$f_{\rm g}=0.5$, $M_{\rm s}=9\times 10^{7}\ {\rm M}_{\odot}$, $R_{\rm s}=1.7$
kpc, $M_{\rm gc}=10^{6}\ {\rm M}_{\odot}$, $R_{\rm gc}=20$ pc, $f_{\rm
rot}=0.3$, inclination angle $\theta_{\rm gc}=30$ degrees, and $R_{\rm
HI}=200$ pc ($D_{\rm HI}=400$ pc). The model parameters for the fiducial model
are summarised in Table LABEL:tab:model_parameters. We also investigate models
with different values of these parameters to discuss the observed diverse
physical properties of GCs. The model parameters for these are summarised in
Appendix B.

Table 1: Description of key physical properties for the fiducial model. Parameters | Values
---|---
Dark matter mass | $3\times 10^{10}{\rm M}_{\odot}$
$c$ parameter | 16
Virial radius of dark matter halo | 33.6 kpc
Stellar disc mass | $9\times 10^{7}\ {\rm M}_{\odot}$
Gaseous disc mass | $9\times 10^{7}\ {\rm M}_{\odot}$
Stellar disc size | 2.4 kpc
Gas disc size | 2.4 kpc
Mass of 1G stellar system | $10^{6}\ {\rm M}_{\odot}$
Size of 1G stellar system | 20pc
Mass of a GC-hosting GMC | $10^{7}\ {\rm M}_{\odot}$
HI hole radius after 1G formation | 200 pc
Mass resolution (gas) | $900\ {\rm M}_{\odot}$
Size resolution (gas) | 20 pc
Mass resolution (GC) | $5\ {\rm M}_{\odot}$
Size resolution (GC) | 0.17 pc
Number of AGB particles per a gas particle | 5
AGB yield | Karakas 2010
AGB feedback effects | Included
Threshold gas density for star formation |  $10^{3}$ cm-3

##  3 Results

In this work we present an overview of our fiducial model; a 47 Tuc analogue.
We analyse the dynamical evolution and accretion history, spatial
distribution, kinematics and chemical properties of the 1G and 2G populations.
Later, we present our global scaling relations of both cluster and galactic
parameters and examine the variety of outputs generated via these simulations.

###  3.1 47 Tuc Analogue

We chose the Galactic GC 47 Tuc as a benchmark for our simulations due to the
detailed and comprehensive observational studies available for this cluster.
We note that several different combinations of parameters resulted in 47 Tuc-
like features (e.g. high 2G mass fraction, evidence of rotation), however, we
summarise only one example here. As 47 Tuc is considered to be a metal rich
cluster, we adjust the parameters of our host galaxy to be representative of a
lower redshift galaxy. We summarise the parameters of the parent galaxies
initial conditions in Table 1. For the GC progenitor, we use a high mass
initial 1G capable of generating the expected He and $\alpha$-element
enrichment patterns (although numerical quantifications of the elemental
composition of the two populations are not considered for our simulation) and
a HI hole size relative to the power emitted from SNe from the 1G. After 370
Myr of evolution, the cluster undergoes significant changes, namely, the
creation of a 2G. Table 2 summarises the initial and final values for our
fiducial model. We acknowledge that there are significant discrepancies
between our stated values and the present-day mass of 47 Tuc. However, the
short duration of our simulation leaves $\approx 10$ Gyrs of evolution between
our masses and that of 47 Tuc. Hence we overestimate both the 1G and 2G masses
to account for any dynamical processes which could result in mass loss over
this period.

####  3.1.1 Dynamical Evolution

Figure 1: The mass of each component (excluding dark matter) in the simulation
as a function of time for the fiducial model. The 1G cluster starts with a set
mass of $10^{6}\rm{M}_{\odot}$ and is embedded within a $10^{5}\rm{M}_{\odot}$
population of disc stars. The HI hole induces a 10 Myr delay in the accretion
of ISM from the galaxy. Following the commencement of gas accretion, there is
a sharp jump in 2G mass associated with both star formation within the 1G and
clumpy accretion. The total GC component represents the sum of the 1G and 2G
masses.

Here we describe the dynamical evolution of a GC (47 Tuc analogue) after its
formation in an exponential gas disc with very clumpy ISM. The simulation runs
in two phases. In the first phase, we evolve a galaxy with an exponential disc
for a period of 100 Myr. This is to create gas density perturbations to select
possible cluster progenitors. After this time period, a high density gas clump
and its surrounding region is replaced with a cluster to approximate the
formation of the 1G. Apart from the selection of a high density region, we do
not perform any further analysis on this phase in the present investigation.
The second phase is the focus of our study and henceforth we use T = 0 to
represent the start of this simulation. We place the 1G in the previous
location of the high density GMC and create a $R_{\rm{HI}}$=200 pc HI hole to
emulate SNe effects of the 1G. Thereafter, the simulation records the
evolution of various parameters within a radius of 30 pc from the centre of
the 1G (1.5$\times R_{\rm{gc}}$). Fig. 1 plots the evolution of 1G, 2G and the
sum of these two components to give the total mass of the GC. We decompose the
gas particles into two separate components; ISM and AGB gas. ISM is
representative of pristine gas originating from the parent galaxy whereas AGB
ejecta is enriched gas released from stars during the AGB phase. Disc stars
originate from the galactic disc but have settled within the cluster’s
potential.

At time T = 0, there is a $10^{6}\rm{M}_{\odot}$ GC progenitor (as stated in
Table. 2) and $10^{5}\rm{M}_{\odot}$ of disc stars within the radius of the
cluster. The 200 pc HI hole causes a $\approx 10$ Myr delay before the GC
starts to interact with gas from the galaxy. Gravitational attraction is the
primary mechanism by which the 1G accumulates the ISM. At a similar time, AGB
stars release low-density gas into the GC. As the 1G stars are the source of
the AGB ejecta, this gas accumulates in the centre of the cluster whereas the
ISM particles distributes itself at a larger radii. Within 50 Myr, the mass of
the 2G is over 50% of its final value. The initial spikes in 2G mass (e.g. at
T=35 Myr, T=45 Myr ) are due to what we have termed ‘clumpy accretion’. Clumpy
accretion occurs when new stars form externally to the GC but fall into its
gravitational potential and mix with the existing population. We observe this
phenomenon in most simulations, so we discuss the impact on the chemistry of
the two populations later in this work. In the first 50 Myr, the mass of the
2G reflects the mass of the ISM, which could illustrate a coupling between the
ISM and the formation of new stars. Later, the introduction of AGB gas becomes
significant and the mixture between pristine and enriched gas creates the
remainder of the 2G stars. Disc stars from the galaxy play a minor role in the
dynamical evolution of the GC; maintaining a relatively constant mass of
$10^{5}\rm{M}_{\odot}$. However, if the gravitational potential of the cluster
trapped these disc stars, it may increase the mass of the GC after 350 Myr by
$\sim 3\%$ based on the GCs 20 pc radius. These disc stars could be a
potential candidate for the older population of 47 Tuc seen though the
photometry of the sub-giant branch (see Milone et al. 2012). The mass of disc
stars within the cluster after the long term dynamical evolution is still an
open question. We intend to investigate this as well as the probability of
natal GMC’s capturing disc stars in future works.

Figure 2: The orbital evolution of the cluster from T = 0 to T = 370 Myr. The
coloured dots are placed at 10 Myr intervals and range from red at T = 0 to
blue after 370 Myr. The gas surface density distribution of the host galaxy at
the start of the simulation has been drawn under the orbital trajectory for
reference.

In Fig. 2, we plot the orbital path of the GC progenitor within its host
galaxy. The X-Y axis defines the spatial extent of the clusters traversal of
the galaxy, with the path taken by the GC shown with a solid black line. A red
star delineates the clusters starting position, and we place a rainbow point
at 10 Myr intervals in order to identify the clusters location at a specific
time. The cluster starts $\approx 85$ pc from the galactic centre and has an
initial velocity of $26\ \rm{kms}^{-1}$. The main accretion events occur from
20 to 50 Myr (i.e. Fig. 1). During this time, the progenitor cluster
experiences a substantial increase in mass and interacts with neighbouring
high density regions. A combination of centrifugal forces and sling-shot
effects from neighbouring clusters forces the GC progenitor into a larger
orbit. The cluster settles into a stable orbit and reaches an apocenter of
$\approx 700$ pc, after 300 Myrs. We plot the initial distribution of gas to
illustrate the various density perturbations as well as the scale of the HI
hole. The high density filaments of gas are the primary source of fuel for
future star formation within the cluster.

Complementing the previous plots, Fig. 3 presents the evolution of new stars
which form during the simulation. The first time step in the top left corner
contains a white circle which denotes the starting location of the 1G. The
next time step occurs after the collapse of the HI hole where the 1G has
established a population of stars from a mixture of AGB ejecta and pristine
gas. We see this population at 71 Myr in the top centre panel, with the white
circle encompassing a small population of new stars which make up the
foundations of the 2G. The remainder of the panels show this population
evolving through the simulation.

The initial distribution of both the GC and the host galaxy greatly influence
the development of the GC. The companion plots to Fig. 3 in Appendix A
illustrate the dynamics of gas particles and disc stars. The first phase of
the simulation generates the turbulent patterns which spawn GMCs observed in
Fig. 15 and allows us to select a gas cloud with mass $>10^{7}\rm{M}_{\odot}$.
At T = 0, high density regions of gas and 2G stars coincides with over-
densities of disc stars. Fig. 1 indicates that the total mass of disc stars
trapped at the commencement of the simulation is $\approx
10^{5}\rm{M}_{\odot}$. Both Fig. 1 and Fig. 16 demonstrate that this mass
remains relatively constant over the next 370 Myr. Our simulation provides
evidence that GMCs which form within a galaxy may have the potential to
capture a population of disc stars before the 1G of a GC is established.
Furthermore, these older disc stars may constitute the metal poor anomalous
population observed in some Galactic GCs (e.g. Terzan 5; Ferraro et al. 2009).
We discuss this further in Section 4.6 and intend for it to be the focus of
forthcoming papers in this series.

In the present simulations, the first phase of evolution restricts gas
particles from transforming into new stars. However, more physical simulations
would allow for high density regions of gas to transform into disc stars. With
this prescription implemented however, the outcome of these simulations were
nearly identical to that of our fiducial model.

Figure 3: The spatial distribution of new stars which formed during the simulation. The X-Y projection spans a 3 kpc range, centred on the parent galaxy. The text in the top left hand corner denotes the time at which that snapshot was taken. In each of the panels, a white circle denotes the location of the 1G progenitor at each time step. After 71 Myr, a number of low mass clusters emerge out of high density gas perturbations and these continue to evolve throughout the duration of the simulation. The corresponding 1G, gas and disc star figures are located in Appendix A.  Table 2: The parameters at the start and end of the 370 Myr evolution of our fiducial model. This model was chosen as a potential candidate for modelling the evolution of the Galactic GC 47 Tuc. The discrepancy between these masses and 47 Tuc’s present day mass is to account for the $\approx 10$ Gyrs of evolution after the creation of the 2G. |  | Initial |  | Final
---|---|---|---|---
$\rm{M}_{1G}<R_{5}$ $(\rm{M}_{\odot})$ |  | 5.07$\times 10^{5}$ |  | 3.95$\times 10^{5}$
$\rm{M}_{1G}<R_{10}$ |  | 8.51$\times 10^{5}$ |  | 6.47$\times 10^{5}$
$\rm{M}_{1G}<R_{20}$ |  | $10^{6}$ |  | 7.37$\times 10^{5}$
$\rm{M}_{2G}<R_{5}$ |  | - |  | 4.31$\times 10^{5}$
$\rm{M}_{2G}<R_{10}$ |  | - |  | 1.10$\times 10^{6}$
$\rm{M}_{2G}<R_{20}$ |  | - |  | 2.19$\times 10^{6}$
$\rm{V}_{x,\ 1G}$ $(\rm{km\ s}^{-1})$ |  | -3.77 |  | -1.65
$\rm{V}_{y,\ 1G}$ |  | 3.77 |  | -1.66
$\rm{V}_{x,\ 1G}$ |  | - |  | 0.03
$\rm{V}_{x,\ 2G}$ |  | - |  | -9.02
$\rm{V}_{y,\ 2G}$ |  | - |  | 6.98
$\rm{V}_{x,\ 2G}$ |  | - |  | -0.31
$\sigma_{x,\ 1G}$ |  | 9.62 |  | 11.04
$\sigma_{y,\ 1G}$ |  | 9.63 |  | 11.07
$\sigma_{z,\ 1G}$ |  | 10.51 |  | 10.84
$\sigma_{x,\ 2G}$ |  | - |  | 6.47
$\sigma_{y,\ 2G}$ |  | - |  | 5.34
$\sigma_{z,\ 2G}$ |  | - |  | 2.61

####  3.1.2 Morphological Properties

Figure 4: The spatial density relation for the 1G, 2G, gas and disc stars
after 370 Myr of the simulation. The 1 kpc box shows the complex structure
surrounding the GC. High-density clusters are visible in the bottom left
corner for the 2G, gas and disc star components.

Fig. 4 shows the final surface density distributions of the 1G, 2G, gas and
disc stars. As expected, the 1G component appears as a single, high-density
point within the simulation box. The 2G component has fragmented into multiple
clumps with a variety of sizes and masses. The density perturbations in the
gas component, which represents both the ISM and AGB gas, mirrors the
distribution of the 1G and 2G particles. The panel depicting the gas component
also illustrates that there is a non-negligible mass of gas remaining within
the cluster. As present-day GCs are known to be incredibly gas-poor, there
must be some mechanism for expelling this gas during the remainder of the
cluster’s evolution. Possible mechanisms are briefly discussed in Section 4.3.
The over density of disc stars seen in the very first time step still surround
the cluster. The 2G population has fragmented into multiple clusters with
comparable sizes and densities to our fiducial cluster. We observe clusters in
the lower left quadrant of the gas and 2G panels undergoing a merging event.
Additionally, several smaller clusters are in the process of being accreted
into larger ones. This points to a hierarchical nature of the accretion
process.

The structure and distribution of the 1G and 2G are well documented in the
literature. A well known criteria of theoretical models of multiple population
formation is that the 2G is more centrally concentrated than the 1G. We
investigate this for our models by looking at the morphologies of the 1G and
2G in a 40 pc box surrounding the GC. Fig. 5 is centred on the 1G’s centre of
mass and illustrates the spatial extent of the 1G and 2G components. The
majority of the 1G’s mass is localised within a 10 pc radius with contours
drawn at $5\times 10^{2}\rm{M}_{\odot}$, $1\times 10^{3}\rm{M}_{\odot}$ and
$5\times 10^{3}\rm{M}_{\odot}$. From the X-Z projection, the 1G has been
slightly elongated due to its angular momentum but is approximately
spherically symmetric. In the X-Y projection of the 2G stars, the central
isophote at $5\times 10^{3}\rm{M}_{\odot}$ occurs at an equivalent radii to
the 1G. However the $5\times 10^{2}\rm{M}_{\odot}$ and $1\times
10^{3}\rm{M}_{\odot}$ isophotes illustrate the complex distribution of 2G
stars. Firstly, the extended structure located on the right side of the figure
is indicative of interactions with ex-situ 2G stars from the environment. The
gravity from the 1G and 2G structures has the potential to tidally attract and
disrupt smaller 2G clusters within the parent galaxy. The spatial extent of
the 2G is much larger than the 1G. One striking difference between the 1G and
2G is the disc like structure of the 2G seen in the X-Z projections. The disc
forms parallel to the plane of the galaxy as a result of its angular momentum.
Mastrobuono-Battisti & Perets (2016) investigated the long-term evolutionary
effects of 2G stellar discs in dense star clusters using N-body simulations
and found that as the 2G disc relaxes, the 1G stellar population flattens the
GC structure to become more elliptical. In our simulations however, the 2G
disc extends far beyond the radius of the 1G and could potentially be tidally
stripped before the relaxation process is complete. We intend to study the
long term evolution of this system in future works to investigate the final
distribution of the 1G and 2G stars within the cluster. The colour map shows
that the central regions of the 1G and 2G have comparable densities. Equal
central densities are observed in some GCs, however Milone et al. (2012)
showed that the central 5 pc had a ratio of $\approx 0.8$ 2G stars. This ratio
decreased to $\approx 0.6$ within 15 pc (conversions from arc minutes to
parsecs used the distance values from Harris 1996(2010 edition) ). Given the
unrelaxed nature of the 2G, future studies will test whether this extended
halo of 2G stars collapses into the centre, thus causing a higher
concentration akin to 47 Tuc.

Figure 5: The X-Y and X-Z distributions of the 1G and 2G populations. Smoothed
contour lines are drawn at $5\times 10^{2}\rm{M}_{\odot}$, $1\times
10^{3}\rm{M}_{\odot}$ and $5\times 10^{3}\rm{M}_{\odot}$ to emphasise the
asymmetry of the 2G. Both the distribution and contours have been smoothed
using a Gaussian filter function.

As the newly formed GC system is in a highly disturbed state, it would be
naive to compare observational radial distributions directly to our results.
However, we include Fig. 6 for completeness. As demonstrated by Fig. 5, the
central concentrations between the 1G and 2G are comparable in mass. Fig. 6
plots the fraction of 2G stars, $f_{2G}$ out to a radius of 20 pc. Within the
central 5 pc, the fraction of 2G stars (i.e. $\frac{M_{2G}}{M_{2G}+M_{1G}}$)
is slightly below $\approx$ 50%. As the density of the 1G decreases, the
extended nature of the 2G results in the outer regions of the GC being
entirely dominated by the 2G. In their current state, these results align with
observations from some clusters which have a centrally concentrated 1G but not
with those from 47 Tuc. Again we highlight that our analysis is based on a
newly born 2G and that its evolution across time may significantly affect its
spatial distribution within the cluster. The results of evolving this system
further in time will be presented in future studies.

One repercussion of Fig. 6 is the prediction that very young GCs should have a
colour gradient. Newly formed GCs should become bluer in colour at increasing
radii due to the higher fraction of younger stars. This appears counter-
intuitive as 2G formation should be predominantly taking place within the
centre of the 1G. We intend to investigate the timescale in which these stars
are either accreted or tidally destroyed by the host galaxy in future studies
of the long term evolution of the cluster. Observational confirmation of any
colour gradients in high redshift clusters may be possible in the future using
high resolution photometric and spectroscopic instruments.

Figure 6: The 3D radial density distributions for the 1G and 2G stars and
corresponding fraction of 2G stars. In the central 5 pc of the cluster, the 1G
and 2G densities are approximately equivalent. However at higher radii, the 2G
becomes the dominant source of mass within the cluster. The shaded region
represents the $1\sigma$ standard deviation.

####  3.1.3 Kinematic Properties

Figure 7: Rotation amplitude as a function of $\theta$ for the fiducial model.
The red and blue shaded regions corresponds to the $1\sigma$ velocity
dispersion. The amplitude of the 2G is far more significant than the amplitude
of the 1G. The GC is rotating around the z-axis and thus we do not expect any
rotation within this component.

Primarily driven by Gaia DR2, kinematical studies of GCs have recently
experienced a growth in popularity. As such, we investigate the kinematics and
internal rotation of our fiducial model. The simulation includes the option to
induce rotation into the 1G population. Additionally, the GC can build up
angular momentum from the accretion of both stellar clumps and gas, as well as
gaining some residual rotation from its orbit around its galaxy. As evident
from Table 2, the process of accreting new material onto the progenitor
cluster drastically alters its kinematics. Milone et al. (2018) used Gaia DR2
to analyse the rotation of the 1G and 2G for 47 Tuc and found that the
rotation of the two populations was in-phase and had similar amplitudes. Fig.
7 plots the amplitude of velocity in the x, y and z directions as a function
of the angle from positive X axis, $\theta$, for our simulated 47 Tuc
analogue. The velocities were obtained by dividing each population into 16
equal-sized bins ranging from $\theta$ = 0 to 360. The average x, y and z
velocity were calculated for each bin and plotted against the corresponding
mean value of $\theta$. The shaded region illustrates the velocity dispersion
for each bin. Next, we apply a sinusoidal curve to each population and compare
the derived amplitudes. The rotation of our two populations is in phase, which
is in agreement with observations of 47 Tuc. However, the amplitude of the 2G
is much larger than that of the 1G. Furthermore, as demonstrated by the shaded
regions denoting the velocity dispersion, the 2G exhibits far more coherent
rotation compared to the 1G, which is at odds with most observations of GCs.
In a study of M80, Kamann et al. (2020) found that the more nitrogen-enriched
population rotates faster than the other two populations which aligns with our
results. Additionally, the aforementioned study by Mastrobuono-Battisti &
Perets (2016) states that their 2G stellar population is characterised by a
lower velocity dispersion and a higher rotational velocity compared with the
primordial older population which is echoed by these results. In our future
works investigating the long term evolution of the cluster, we intend to
explore whether the differences between the rotational velocity of the 1G and
2G diminishes as to align with current observations of Galactic GCs. Using a
suite of N-body simulations, Hénault-Brunet et al. (2015) determined that
different pollution scenarios could result in distinct kinematic properties
which could be used to distinguish between various scenarios. Similarly, we
intend to analyse the kinematic imprint of the different rotation amplitudes
between each population for our scenario, the outcome of which could either
strengthen or falsify our model.

Figure 8: The velocity map along the line of sight for the 2G. Contours are
placed at densities of $1\times 10^{3}\rm{M}_{\odot}$, $5\times
10^{3}\rm{M}_{\odot}$ and $1\times 10^{4}\rm{M}_{\odot}$ to show the spatial
extent of the stars. The rotation map has undergone Gaussian smoothing however
the contours are as-is to show additional structure within the plane of the 2G
disc.

Turning our attention to the rotation of the 2G only, Fig. 8 shows the
velocity map along the line of sight. A Gaussian filter has been applied to
smooth over any outlying escaping 2G star particles. The contours serve as a
reference for the mass isophotes. The 2G disc is clearly rotating along the
line of site with a maximum velocity in the $y$ direction of $\approx
8\rm{kms}^{-1}$. The high velocity dispersion and low rotation of the 1G cause
any signs of rotation to be far less obvious and thus this plot is omitted
from the present study.

####  3.1.4 Chemical Enrichment

Studying chemical compositions in GCs is central to their characterisation.
Clusters which are homogeneous in heavy element abundances are known as Type I
GCs whereas Type II GCs are those with more than two dominant populations. Our
simulations can reproduce both types of GCs through altering the time, mass
and starting location of any accreted clumps.

Figure 9: The helium abundance for the 2G for three different $Y_{AGB}$
values. The shaded regions represent the $1\sigma$ standard deviation. The
transition to the helium abundance of pristine gas ($Y=0.25$) corresponds to
the radial extent of the 1G.

One of the main strengths of the AGB scenario is its ability to explain the
Mg-Al anti-correlation and other elemental patterns, as well as being able to
deposit a significant mass of Helium (He) into the GC. We estimate the He
enrichment of the 2G via the formula:

$Y_{2G}=\frac{Y_{AGB}\times M_{AGB}+Y_{ISM}\times M_{ISM}}{M_{AGB}+M_{ISM}}$
(5)

where $Y_{AGB}$ is the typical He abundance for AGB ejecta and $Y_{ISM}$ is
the He abundance for pristine gas in the ISM. $M_{AGB}$ and $M_{ISM}$ are the
masses of new stars formed from AGB ejecta and pristine ISM respectively. Fig.
9 shows the He enrichment as a function of radius with the shaded regions
corresponding to the $1\sigma$ standard deviation. We assume three different
He abundances in AGB gas; $Y_{AGB}=$ 0.3, 0.34 and 0.38, and pristine ISM to
be $Y_{ISM}=0.25$. Stars in the central 9 pc are He enriched, with
$Y_{2G}>Y_{ISM}$ for all He abundances. This is expected as this region has a
high density of 1G stars (Fig 6) which are the primary source of AGB ejecta.
After 9 pc, there is a steep decline in the He enrichment and at radii greater
than 12 pc, the He abundance drops to approximately that of $Y_{ISM}$. This
calculation assumes that $Y_{ISM}$ is constant throughout the galaxy when more
realistically, the GC could gravitationally interact with inwardly or
outwardly moving gas in the disk. This could allow the GC to accrete gas with
various $Y_{ISM}$, generating a larger $\Delta Y$ between the 1G and 2G.

An early study of 47 Tuc by Briley (1997) suggested that a non-negligible
spread in the He abundance may explain the morphology of the horizontal
branch. Later, Anderson et al. (2009) required a $\Delta Y=0.026$ in order to
explain the colour spread in the main sequence, provided that He is the main
source of broadening. Today, the observed He abundance of 47 Tuc is given to
be $Y_{1G}=0.256$ and $Y_{2G}=0.276$ (Fu et al., 2018). Super AGB stars (e.g.,
Siess 2010; Ventura & D’Antona 2011; Doherty et al. 2014) may be capable of
producing a He abundance of $0.38$, however this is still not enough to
generate the $\Delta Y$ observed for 47 Tuc. Approaches to this problem
include increasing the He concentration of the pristine gas to match the 1G or
decreasing the dilution of AGB gas by limiting ISM accretion. Additionally,
due to the uncertainties around AGB yields, future revisions of this model
with updated AGB prescriptions may be able to rectify this difference.

Figure 10: The radial distribution of in-situ and ex-situ stars within 1 pc
bins. Stars formed in-situ to the 1G make up the majority of the central
region of the 2G. Stars which formed externally to the 1G make up the dominant
population outside the 10 pc radius.

The occurrence of in-situ and ex-situ stars has yet do be discussed in the
context of GC evolution. Clumpy accretion within proto-GCs naturally explains
Type II GCs, but the chemical consistency of Type I GCs may pose some issues
for our scenario. Furthermore, there is no guarantee that these new stars are
genuine members of the 2G population and may even have chemistry
representative of the 1G. Given the degree of He enhancement of the 2G within
the central region of the GC, we would expect that a significant proportion of
these stars were formed within the GC environment. We confirm this hypothesis
in Fig. 10, which plots the mass of in-situ and ex-situ stars at increasing
radii. Within the central ten pc which we know is dominated by the 1G, the
majority of the stars formed in-situ with the GC. With increasing radius, the
dominant population of 2G stars becomes those which were formed outside the
1G, and were assimilated into the GC at a later time during a clumpy accretion
event.

Figure 11: The probability density plot of the original location of new stars.
The in-situ and ex-situ fractions have been scaled according to their total
contribution to the total 2G. The original locations of each type of new star
peaks at the same location and thus should have similar metallicities.

In addition to variations in He, Lithium (Li) abundances also poses a
challenge to many MSP formation scenarios (for a recent discussion, see
Gratton et al. 2019). Forged during the Big Bang nucleosynthesis, Li may be
vital in determining the stellar source of GC polluters. The AGB scenario has
been heralded as a promising solution to the problem of Li variations between
populations due to AGB stars ability to manufacture Li through the Cameron
Fowler mechanism (Cameron & Fowler, 1971). However, the abundance of Li in
pristine material accreted onto the GC has been suggested to influence this
result (D’Antona et al., 2012). Here we do not consider the influence of
lithium production and dilution effects due to the vast uncertainties from AGB
yields; however, our model’s ability to generate a diverse range of 2G stars
enriched by various quantities of AGB gas works in our favour.

The metallicity gradient of galaxies is well documented in the literature
(e.g., Searle 1971; Magrini et al. 2016; Bresolin 2019). We assume that if the
accreted clump originated from a similar location within the GC, it would have
similar metalicities to the accreted ISM. Fig. 11 shows the distribution of
the original locations of all 2G star progenitors which finished within a
radius of 20 pc from the centre of the 1G in the final time step. The navy ex-
situ curve and gold in-situ curves are scaled according to their total
contribution to the 2G. With respect to the centre of the parent galaxy, the
maximum value of the probability density distribution demonstrates that the
majority of 2G material originated from $\approx$ 500 pc from the galaxies
centre. This is reasonable given the GCs initial offset from the centre and
radius of the HI hole, as seen in Fig. 2. This figure shows that there is a
significant spread in the initial positions of the gas accreted by the GC,
with the maximum value of an ex-situ particle which made its way into the GC
at 1750 pc. This introduces a potential problem for our scenario; gas
particles from such a large radii may introduce a non-negligible [Fe/H] spread
which is not detected in Galactic GCs. However, as we are analysing the
starting positions of all new stars within a 20pc radius, Fig. 9 and 10
demonstrate that it is questionable whether anything at radii larger than 10
pc is a genuine 2G star. Future investigations into the long term evolution
will again affect this distribution and whether this spread in initial
locations of 2G stars is a major issue for this model.

###  3.2 Scaling Relations

The results given in Milone et al. (2020) demonstrate the variety of $f_{2G}$
vs GC mass relations which exist for GCs and how their parent galaxy can
influence the gradient of the slope. Here we present similar scaling
relations, derived from our simulations. The results should be taken with
caution as the dynamical effects such as stripping and SNe have the potential
to dramatically influence the correlation between the 2G and total mass of a
cluster. Here we summarise the most notable correlations using a data set of
simulations which is representative of Galactic GCs. 34 clusters meet the
criteria which include having a realistically sized HI hole, thus allowing gas
accretion to occur within 370 Myr. For the present simulations, this
translated into a radius of no more than 300 pc. We also require that the 1G
is not destroyed in some way by the parent galaxy. This results in setting a
minimum size threshold on possible progenitors. Several simulations ran for
only 150 Myr but have been excluded as the proto-GCs are usually in a far more
turbulent state, and for some models, have not completely finished their
accretion process. We also provide a sample of AGB only simulations which we
achieve by modifying the HI hole to be larger than 500 pc, artificially
inhibiting any gas accretion, which we use as a comparison for our accreting
GC models. Finally we present a sample of GCs formed in younger galaxies or
low surface brightness galaxies (LSB), both of which are gas poor when
compared to our fiducial model’s galaxy. The parameters for each of these
models are available in Tables 4, 5 and 6 in the Appendix.

####  3.2.1 The $f_{enriched}$ vs $M_{GC}$ relation

Figure 12: $f_{1G}$ as a function of GC mass within a 5 pc, 10 pc and 20 pc
radius. The line of best fit, shown by a solid grey line, is calculated using
only the Galactic GC data. The equations of these lines can be found in Table.
3. The dashed line in the bottom panel is representative of the initial masses
predicted by Baumgardt et al. 2019, taken from Fig. 6 in Gratton et al. 2019.

The relationship between the mass of the 2G and the total mass of the cluster
is one of the most notable results arising from the last decade of GC
research. From our simulations we recover a correlation which agrees with the
results from observations. To be consistent with prior works in the
literature, in Fig. 12 we plot this relationship relative to the fraction of
1G stars within the cluster, $f_{1G}$. We expect to find small clusters
dominated by 1G stars residing in the top left, and more massive 2G dominated
clusters in the bottom right. The values represent the initial masses of
possible Galactic GC progenitors immediately after the 2G has been formed. For
the Galactic GC models, the 1G half mass radius was approximately at 5 pc. We
do not consider the radius of the 2G as the stars tend to be in a turbulent
state. For all stars within a 5 pc cut off, there is a moderate anti-
correlation between the total mass of the cluster and fraction of 1G stars
(see Table 3). A steeper gradient may exist for clusters within the mass range
of $10^{6.5}\rm{M}_{\odot}$ to $10^{5.5}\rm{M}_{\odot}$, however outlier
points on the lower mass end mask this result. Increasing the cut off radius
to 10 pc, which in most cases contains the majority of 1G mass (see Fig. 6),
this alternative steeper gradient still exists but the correlation between the
mass of the cluster and the fraction of 1G stars improves. Finally, we see our
strongest correlation at a radius of 20 pc. A useful comparison for these
values is initial masses of the Galactic GCs presented in Baumgardt et al.
(2019). The data clusters at a slightly higher GC mass than the results
predicted by Baumgardt et al. (2019); however, we expect this is more of a
selection effect of our study primarily focusing on large, 47 Tuc-like
analogues. Our fiducial model discussed in the previous section lies on top of
the predicted curve. Previous investigations into GC formation scenarios have
emphasised that pristine gas accretion is necessary in explaining the O-Na
anticorrelation. Our study adds to this argument that gas accretion is
required in order to reproduce the relationship between the initial mass of
the cluster and the fraction of 2G stars. However, this result relies on the
assumption that all 2G stars within a 20 pc radius are accreted by the GC. Fig
5 and 6 illustrates that at least for the fiducial model, the outer regions of
2G stars may be preferentially stripped from the cluster during further
evolution. As this is a determining factor in the success of this model, we
intend to study the long term evolution of the cluster in more detail in
future works. Our model emulates Magellanic Cloud GCs or GCs which originated
in younger galaxies by simulating Low Surface Brightness (LSB) galaxies. The
clusters from these simulations always appear more 1G dominated than their
Galactic counterparts and lie towards the right of the line of best fit. This
agrees with the results presented in Milone et al. (2020) that Magellanic
cloud clusters have a lower fraction of 2G stars than Galactic GCs.

Fig. 12 also illustrates the lower limit on the mass of the 1G in order for a
realistic 2G to form. M21 in Table 4 represents the lowest mass GC that was
able to retain a 2G population. As the size of the GC decreases, it starts to
become more susceptible to the gravitational influence of larger clusters
which formed within the Galaxy. However, when in a LSB galaxy, the lower limit
for 2G formation differs as there are fewer massive clusters which disturb the
gas around the 1G. This can be seen for M44 in Table 6, particularly within
the central 5 pc of the cluster. Here we would predict a much lower 2G
fraction for that mass, however, altering the parameters of the host galaxy
allows this model to account for a wide range of observations.

Figure 13: The mass of the 2G against the total mass of the cluster within 20
pc. The Galactic GC points are tightly correlated. AGB only models do not
follow this trend. The coefficients extracted from fitting lines to these data
sets can be found in Table 3.

We exclude Bondi accretion (Bondi, 1952) as the main mechanism by which gas is
accreted onto the GC. To confirm this, we plot the mass of the 2G against the
mass of our Galactic GCs in Fig. 13. The coefficients listed in Table 3
demonstrates a strong correlation between the 2G and total GC mass and that
the corresponding gradient is approximately 1.5. Bondi accreation assumes that
the amount of accreted material is proportional to the square of the stellar
mass. Our derived slope is not as steep as that predicted by spherically
symmetrical accretion, and we assume that gas accretion is due to the
gravitational potential of the 1G. In comparison to this, Table 3 also
includes the results from fitting our AGB only models. Although there is a
much smaller sample size, the line of best fit has a gradient of 2.55. We
expect this to be different to our Galactic GC models as no gas accretion is
involved.

####  3.2.2 The rotation amplitude vs $M_{GC}$ relation

One unexpected correlation found during the development of this scenario was
the connection between the final mass of the cluster and the rotation
amplitude of 2G stars. By extracting the amplitude from both populations (see
Fig. 7), Fig 14 compiles the results from our physical simulation samples to
illustrate that the rotation of the 2G is predicted to be higher for more
massive clusters. Table 3 lists a moderate positive correlation for $V_{x}$
and $V_{y}$ vs $\rm{log_{10}(M}_{GC})$ parameter pairs. No correlation exists
for the $V_{z}$ \- $\rm{log_{10}(M}_{GC})$ parameter pair as this represents
the axis of rotation of both the cluster and the parent galaxy. Similarly, the
difference between the maximum rotation amplitude of the 1G and 2G exhibit a
moderate positive correlation. We expect this result to be invariant to the
host galaxy parameters as this trend is also evident in our LSB models. This
finding conflicts with current observations of Galactic GCs which, for the
most part, show little to no difference between the rotation amplitudes of the
two populations (e.g. see Cordoni et al. 2020). Furthermore, this relation
implies that larger clusters will exhibit signs of internal rotation. Current
studies of internal kinematics using Gaia DR2 have shown this not to be the
case, although future data releases with more precise measurements may alter
this view. It is unclear how the dynamical evolution of our simulated GCs will
progress given the different rotation characteristics. We hope that the long-
term evolution can significantly weaken the rotation amplitudes of the two
populations and we intend on investigating this question in future studies

Another kinematic prediction is that even when the initial internal axis of
rotation ($\theta_{\rm{gc}}$) for the 1G is altered, the accretion processes
forces the rotation axis of the proto-GC to align with the parent galaxy. For
example, M27 (see Table 4) initially starts rotating at an angle of
$60^{\circ}$ with respect to its host galaxy and after the 2G formation
process, rotates in line with the galaxy. The corresponding AGB only
simulation, M40, retains this initial internal rotation for both the 1G and
2G. In most cases, the 2G in AGB only models share a similar phase and
amplitude to the 1G. This is an important conclusion for our scenario as it
suggests that the host galaxy leaves a kinematic imprint on the cluster. We
discuss this point further in Section 4.5.

Figure 14: Final mass of the GC as a function of the rotation amplitude in the X, Y and Z directions. The coefficients for the lines of best fit are found in Table. 3.  Table 3: Correlations seen between different sets of parameters for the Galactic GC sample. Correlations are determined by fitting the line $y=ax+b$ to the corresponding data sets. Unless specified, the clusters were set to an arbitrary 20 pc cut off. (∗ corresponds to our AGB only model.) Parameter Pair | a | b | r | p-value
---|---|---|---|---
$\rm{f}_{1G}$ \- $\rm{log}_{10}(\rm{M}_{GC})$ < 5 pc | -0.504 | 3.418 | -0.729 | 6.90e-07
$\rm{f}_{1G}$ \- $\rm{log}_{10}(\rm{M}_{GC})$ < 10 pc | -0.449 | 3.137 | -0.768 | 7.02e-08
$\rm{f}_{1G}$ \- $\rm{log}_{10}(\rm{M}_{GC})$ | -0.403 | 2.851 | -0.868 | 1.47e-11
$\rm{log}_{10}(\rm{M}_{2G})$ \- $\rm{log}_{10}(\rm{M}_{GC})$ | 1.477 | -3.211 | 0.987 | 4.08e-27
$\rm{log}_{10}(\rm{M}_{2G})$ \- $\rm{log}_{10}(\rm{M}_{GC})^{*}$ | 2.553 | -11.442 | 0.984 | 1.05e-05
$V_{x,2G}$ \- $\rm{log}_{10}(\rm{M}_{GC})$ | 7.544 | -39.236 | 0.798 | 9.50e-09
$V_{y,2G}$ \- $\rm{log}_{10}(\rm{M}_{GC})$ | 7.834 | -41.342 | 0.713 | 1.55e-06
$V_{z,2G}$ \- $\rm{log}_{10}(\rm{M}_{GC})$ | 0.271 | -0.966 | 0.229 | 1.87e-01
$V_{x,2G-1G}$ \- $\rm{log}_{10}(\rm{M}_{GC})$ | 9.001 | -45.941 | 0.800 | 8.39e-09
$V_{y,2G-1G}$ \- $\rm{log}_{10}(\rm{M}_{GC})$ | 7.996 | -42.379 | 0.713 | 1.56e-06
$V_{z,2G-1G}$ \- $\rm{log}_{10}(\rm{M}_{GC})$ | 0.433 | -2.003 | 0.338 | 4.73e-02

##  4 Discussion

###  4.1 The origin of the $\rm{M}_{2G}$ vs $\rm{M}_{GC}$ relation and the
link to LSB galaxies

The results from Fig. 12 and 13 illustrate how the properties of the parent
galaxy can influence the fraction of 2G stars within the GC. The gas fraction
of the galaxy is the primary factor which mediates the fraction of 2G stars.
We also find a lower limit on the initial mass of a cluster in which we can
establish a 2G population. Clusters which are born closer to the galactic
centre experience more clumpy accretion events due to their densely populated
environment. This results in larger 2G populations for massive clusters,
however there is a higher chance that a cluster may be tidally destroyed by
its siblings. Conversely, clusters which do not pass through the galactic
centre and born at larger radii have a high fraction of in-situ star formation
within the 2G population and are less likely to be tidally disrupted. This
enables low mass GCs to survive and generate a feasible 2G of stars. The
metallicity of the galaxy has minimal consequences on the 2G fraction or
morphology of the resulting cluster, however it does impact the He enrichment,
with higher metallicity galaxies generating slightly more He-rich 2Gs.
Modelling the chemistry of the clusters is difficult with the present
simulations and thus we do not dwell on these results.

As we suggest that the gas density of the host galaxy is one of the primary
drivers of the fraction of 2G stars, we expect that GC formation is not
necessarily a product of the early universe. Observations have demonstrated
that high redshift galaxies are significantly more gas-rich compared to z = 0
galaxies (e.g. Daddi et al. 2010). Young stellar clusters may also be able to
develop GC like qualities provided that they are exposed to extremely high gas
fractions as were common in the early universe. GCs from LMC-like and other
low surface brightness galaxies tend to have a lower fraction of 2G stars when
compared to the Galactic counterparts. As found by Milone et al. (2020), LMC
clusters are more 1G dominated and do not lie on the line of best fit for
Galactic GCs and we see this same phenomenon occurring in our models. A lower
gas density can have the effect of lowering the minimum mass that is needed in
order to create a 2G as illustrated by M44 from Table 6.

###  4.2 The chemical implications of clumpy accretion

The phenomenon of clumpy accretion allows our model to explain several
observable parameters, including the expected mass of 2G stars and a
connection to the parent galaxy’s parameters. However, it introduces several
unknown factors concerning the chemical composition of the cluster. Given the
vast range of radii that the accreted gas particles originated from in Fig.
11, this accretion process may induce unwanted chemical abundance spreads.
Observationally, it is challenging to estimate the [Fe/H] gradient for high-z
GC-hosting dwarf galaxies. Taking the LMC as a proxy, the [Fe/H] gradient is
shallow (-0.01 dex/kpc for old stars, Fig.2 from Feast et al. 2010). Thus if a
GC forms from gas spreading over a 1 kpc region, then Fe/H spread is only 0.01
dex for the 2G alone. An analysis of globular clusters from the APOGEE survey
by Mészáros et al. (2019) stated that the true iron spread for their GCs was
on average 0.068 dex across their 30 clusters. For our fiducial model, gas
accreted to form the 2G spans a radius of 1.75 kpc and depending on the slope
of the [Fe/H] gradient employed for our galaxy, the final spread in
metallicity could be within the range set by observational results.

Clumpy accretion comes with many more obstacles than inducing a spread in
[Fe/H]. By delaying the accretion of 2G clumps, this scenario can naturally
explain Type II clusters such as NGC 2808. However, if clumpy accretion is to
occur for Type I GCs, it would need to take place in the very early stages of
the 2G formation process. Even still, this may induce changes in the chemical
makeup of the 2G. The merging of GCs is discussed in van den Bergh (1996) as a
mechanism for producing composite color-magnitude diagrams. Mergers may have
occurred in dwarf spheroidal galaxies with low internal velocity dispersion.
At T = 0 of the fiducial model, the 1D velocity dispersion of disc stars is
$19kms^{-1}$. This low dispersion may relate to the incidence of clumpy
accretion.

High precision measurements of Galactic GCs have uncovered slight variations
in populations within some regions of a cluster’s colour-magnitude diagram
(CMD). For example, Di Criscienzo et al. (2010) noted that 47 Tuc is comprised
of three different subpopulations; a first-generation consisting of 30% of the
stars in the cluster, a SGI (equivalent to a 2GI in our nomenclature) which
makes up 60% of stars belonging to the second generation and lastly a SGII
(2GII nomenclature) which makes up 10% of 2G stars and is found within the
faint subgiant branch in the CMD. Our model can rationalise the existence of
these 2G subpopulations by attributing them to slight differences in the
location these populations were formed. More enriched populations could
originate within the centre of the GCs where the gas is heavily enriched with
He from AGB ejecta, whereas less He rich stars could form in similar
conditions within the galaxy and then be later accreted by the cluster.
Another example of this phenomenon documented by Milone et al. (2020) is found
in the Small Magellanic Cloud (SMC) cluster NGC121 which also has two enriched
components. Clumpy accretion could still be possible for low surface
brightness galaxies such as the SMC and thus our scenario could potentially
explain this observation. This is all, of course, speculative and requires
both even higher resolution simulations and more precise observations with
which we can compare our results.

###  4.3 2G concentration and long term evolution

Fig. 6 illustrates that the outskirts of the cluster are entirely 2G
dominated. This issue plagues the majority simulations which showed external
gas accretion. In its current state, the distribution of 2G stars around the
cluster poses issues for this scenario as observations show that 2G stars are
generally more centrally concentrated (Lardo et al., 2011). The 2G stars form
in a disc, parallel to the plane of the galaxy rather than showing spherical
symmetry like the 1G. To reproduce observations, the cluster must avoid having
its 2G stars being stripped by the galaxy and redistribute them towards the
centre. However, the Type II clusters M15 and M80 both show a central
concentration of 1G stars (Larsen et al. 2015; Dalessandro et al. 2018) and
thus clumpy accretion may be the cause of these observations. The cores of GCs
with a higher central concentration of 2G stars tend to be more ex-situ star
dominated. As ex-situ stars tend to be less He rich compared to their AGB gas
enriched in-situ counterparts, this would influence any chemical gradients
within the cluster. On the assumption that all 2G stars within a 20 pc radius
of the 1G mass are retained by the GC and that the 2G stars become more
centrally concentrated as the cluster relaxes, this would allow this scenario
to solve the mass budget problem. This result agrees with a study of pristine
gas accretion by Calura et al. (2019), who determined that the degree of 1G
mass loss required to produce observed ratios of 1G to 2G in GCs is less
extreme than what the mass budget problem suggests.

Our fiducial model overestimate the mass of GC progenitor to allow for some
degree of mass loss over its life time. Additionally, there will be a further
reduction of 2G mass due to future SNe events. However, these events will be
necessary to remove the non-negligible amount of gas which remains within the
proto-cluster. For our fiducial model, the simulation ends with $\approx
6\times 10^{4}\rm{M}_{\odot}$ of gas within 20 pc of the centre of mass. This
gas must be expelled in order to explain the lack of neutral ISM observed in
Galactic GCs. Simulations by Chantereau et al. (2020) found that both ram
pressure stripping and ionisation is mandatory to explain the small amount of
ionized gas in the core of GCs. We intend to investigate these channels of gas
expulsion in future simulations of our Galactic GCs.

###  4.4 The physicality of a large HI hole

An evacuated region of gas is placed around the progenitor cluster to account
for SNe effects from the 1G on the surrounding ISM. This HI hole acts to
suppress the early accretion of gas onto the cluster. We admit that this is a
gross simplification, however, instances of HI holes appearing in dwarf
galaxies have been recorded in the literature. Warren et al. (2011) discussed
possible origins of kpc-scale holes in the atomic hydrogen distributions of
some nearby dwarf irregular galaxies. Using radial analysis, they calculate
the HI hole size for DDO 181 to be $\approx 740$ pc, Holmberg I to be $\approx
850$ pc and M81 Dwarf A to be $\approx 745$ pc as three examples. Each of
these radii fit comfortably in the range we set for our cluster progenitors.
However, different methods can result in vastly different predictions of hole
size. Using the THINGS (Walter et al., 2008) survey, Bagetakos et al. (2011)
identified over 1000 HI holes in their sample, and as an example, found six
holes in Holmberg I, ranging in size from 190 pc to 740 pc. Vorobyov et al.
(2004) used numerical simulations to investigate the energy required to create
such a hole in Holmberg I. Through their modelling they deduced that a hole of
this magnitude required an energy equivalent to $300\pm 50$ Type II
supernovae. The energy our 1G is capable of expelling will be tightly linked
to its IMF. A back of the envelop calculation assuming a 1G population with a
Salpeter IMF of $\alpha=2.35$ and mass range of 0.1 - 50 $\rm{M}_{\odot}$
would require $7\times 10^{-3}\rm{M}_{\odot}$ SNII per $\rm{M}_{\odot}$. We
assume $4.2\times 10^{4}\rm{M}_{\odot}$ for HI hole formation and thus with
our $10^{6}\rm{M}_{\odot}$ fiducial cluster, a 200 pc hole could be assumed to
be a conservative estimate. Although the radius of the HI hole ($R_{\rm{HI}}$)
is treated as a free parameter in our model, it is strongly influenced by
feedback effects of the 1G. This in turn impacts the duration of the holes
collapse and thus the subsequent gas accretion onto the GC. We intend to
discuss this further during investigations into 1G formation and HI hole
creation in order to rationalise our present assumptions.

Models which do not include HI holes overestimate the total mass of the
cluster. Secondly, there is no delay in star formation and therefore, no
distinction between 1G and 2G stars. Gas that would have been consumed during
1G formation is readily available to the cluster to immediately start forming
the 2G. There should be some physical mechanism which prevents this from
occurring and thus some form of HI hole is necessary in reproducing
observations.

###  4.5 Spin axis alignment between the GC spin axes and the parent dwarf

As discussed in Section 3.1.3, the accretion process causes the internal
rotation axis of the GC to align with that of the parent galaxy, independent
of the GCs host environment. This prediction has implications for the
identification of accreted clusters and can potentially aid in reconstructing
the accretion history of the Milky Way (MW). Taking the ’Sausage’ GCs
identified by Myeong et al. (2018) as an example, our scenario would predict
that the rotation axis of the integrated orbits of these clusters would align
with the internal rotation of the Gaia-Enceladus dwarf before it was destroyed
by the MW. Additionally, this method could contribute to the authentication of
any potential Sausage clusters. We note that this may only be used for
rotating clusters. NGC 2808, NGC 7089 (M2) and NGC 1904 (M79) are three Gaia
Sausage clusters as identified by Myeong et al. (2018) and have evidence of
rotation as found by Sollima et al. (2019). By integrating the orbits of these
clusters and analysing their axis of rotation at the time they were accreted,
they could be used to test our hypothesis. This is reliant on the long term
evolution not erasing the kinematic signature from the host.

###  4.6 Disc stars as an additional metal poor population

As evident in Fig. 1, the simulation starts with $10^{5}\rm{M}_{\odot}$ of
disc stars surrounding the 1G. The concentration of these stars is visible in
the left most panel of Fig. 4 and as the system evolves this mass remains
relatively constant. Fig. 16 demonstrates that the clustering of disc stars
also occurs for several of the fiducial model’s sibling clusters. This
phenomenon is found across the whole data set, and is not dependent on any
free parameters. From this we conclude that GCs should be capable of retaining
a small population of disc stars from their parent galaxy. The $\sigma_{1D}$
velocity dispersion of the disc stars within a 50 pc radius of the cluster is
$15kms^{-1}$. Our fiducial GC with a mass of $~{}3\times 10^{6}\rm{M}_{\odot}$
and radius of 20 pc could potentially capture these disc stars as the internal
dispersion of the cluster is $\sqrt{\frac{GM_{\odot}}{R}}\approx 25kms^{-1}$.
This may explain the origin of the 3rd metal poor population seen on the
subgiant branch in 47 Tuc (Milone et al. 2012), or the small 3rd generation
seen in Terzan 5 (Ferraro et al. 2009). These captured disc stars are a
feature in all simulations where the 1G progenitor was born in a central
region of the galaxy. Our simulations predict that most Galactic GCs should
comprise of a metal-poor component with the percentage of total mass varying
depending on the dynamical evolution of the cluster. Future works intend on
taking a much deeper look into the formation processes of how a 1G may capture
these stars in the first place and how these stars may evolve after the
cluster has been formed.

##  5 Conclusions

We have presented the preliminary results from our original, hydrodynamical
simulations which builds upon the theoretical framework first proposed by
D’Ercole et al. (2008). In our model, we have first identified GC progenitor
gas clouds with masses of $10^{7}\rm{M}_{\odot}$. Next, these clouds are
converted into compact stellar systems associated with the 1G. Around this
system, we included a HI hole produced by the multiple 1G SNe to demonstrate
how SNe can influence the gas accretion processes and subsequent star
formation. The gas accretion process proceeds after the collapse of this hole
and ISM from the galaxy mixes with AGB ejecta from 1G stars. This resulted in
the birth of a secondary population of stars known as the 2G.

The foremost result from our study is that we have been successful in
producing a correlation between the mass of the cluster and the fraction of 1G
stars. Additionally, the gas fraction of the galaxy is the dominant parameter
in mediating this relationship with high gas fraction galaxies producing GCs
which are more 2G dominated. The gas fraction also sets the lower limit for
which no 2G is observed due to gas poor galaxies within the simulation
experiencing less turbulence which can lead to the disruption of the GC. Prior
to the formation of the 1G, we find that the GMC has the potential to capture
a population of disc stars. We predict that these stars are the missing minor,
metal-poor component that has been observed in some Galactic GCs.
Additionally, we expect that a very small population should be present in the
majority of GCs. We intend to examine this prediction further in our next
paper.

We observe GC progenitors interacting and accreting one another through a
process we have termed clumpy accretion. Clumpy accretion is a hierarchical
process and can occur at any point during the simulation. However, we suggest
that bona fide Type I clusters can only undergo clumpy accretion within the
first 50 Myr as to minimise the risk of inducing a [Fe/H] spread. We argue
that this form of accretion may be beneficial in producing smaller
subpopulations of 2G stars observed in Galactic GCs. Type II clusters on the
other hand, may experience clumpy accretion at any point within the 370 Myr
window of the simulation as this would give rise to distinct populations
characterised by multiple [Fe/H] values. Our scenario suggests that Type I and
Type II GCs share the same origin.

The kinematic repercussion of our scenario is that the 2G of a GC exhibits
highly coherent rotation in a disc-like structure, parallel to the plane of
the parent galaxy. As the long term dynamical effects of this rotation is
unclear, we intend to perform follow up simulations on our cluster progenitors
in future studies. If the GC retains information about the plane of the host
galaxies rotation, this information could be used to inform cosmological
simulations about the accretion history of the Milky Way. Provided this
kinematic imprint of the parent galaxy remains within the cluster, we predict
that by backwards integrating the motion of confirmed ex-situ GCs and
analysing their axis of rotation, it may be possible to infer the plane of the
dwarf galaxy prior to its accretion.

The cornerstone of this model is its ability to reproduce the observed
correlation between the GC’s mass and its fraction of enriched stars. However,
this result relies on the accretion of all 2G stars within a 20 pc radius of
the GC. Further long-term dynamical investigations are required to determine
the outcome of these stars as this could lead to the falsification of our
results. In its current state, this scenario is unable to explain all
observational criteria for Galactic GCs. Nevertheless, it still holds
scientific merit as being the first of its kind to simultaneously model the
dynamics of MSP formation within a GC and the entirety of its parent galaxy.
We intend to continue to improve and develop this model to more closely align
with the observed characteristics of GCs. Gaia DR3 will shed more light on the
internal kinematics of GCs while spectroscopy from instruments such as
VLT/MAVIS and photometry from the James Webb Space Telescope will better
quantify the chemistry of the populations. Further developments in AGB yields
will allow us to refine the degree of mixing with pristine gas to constrain
our parameters further. We have paved the way for self-consistent simulations
of GC formation scenarios and look forward to comparing our results with
future groups who undertake similar investigations.

## Acknowledgements

The authors thank the anonymous referee for their kind and thoughtful
comments. Their suggestions helped to identify ambiguities and improved the
clarity of this work. Numerical simulations were run on the Pleiades and
OzSTAR GPU clusters kindly made available to us through the International
Center for Radio Astronomy Research at The University of Western Australia and
the Centre for Astrophysics and Supercomputing at the Swinburne University of
Technology. This study made use of the PYTHON packages NUMPY (Van Der Walt et
al., 2011), SCIPY (Virtanen et al., 2020), MATPLOTLIB (Hunter, 2007), PANDAS
(McKinney, 2010) and IPYTHON (Pérez & Granger, 2007).

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Anderson et al. (2009) Anderson J., Piotto G., King I. R., Bedin L. R., Guhathakurta P., 2009, ApJ, 697, L58
  * Armstrong et al. (2018) Armstrong B., For B. Q., Bekki K., 2018, MNRAS, 481, 3651
  * Bagetakos et al. (2011) Bagetakos I., Brinks E., Walter F., de Blok W. J. G., Usero A., Leroy A. K., Rich J. W., Kennicutt R. C. J., 2011, AJ, 141, 23
  * Bailin (2019) Bailin J., 2019, ApJS, 245, 5
  * Barbuy et al. (2009) Barbuy B., Zoccali M., Ortolani S., Hill V., Minniti D., Bica E., Renzini A., Gómez A., 2009, Astronomy & Astrophysics, 507, 405–415
  * Bastian & Lardo (2018) Bastian N., Lardo C., 2018, Annual Review of Astronomy and Astrophysics, 56, 83
  * Baumgardt et al. (2019) Baumgardt H., Hilker M., Sollima A., Bellini A., 2019, MNRAS, 482, 5138
  * Bedin et al. (2004) Bedin L. R., Piotto G., Anderson J., Cassisi S., King I. R., Momany Y., Carraro G., 2004, ApJ, 605, L125
  * Bekki (2011) Bekki K., 2011, MNRAS, 412, 2241
  * Bekki (2013) Bekki K., 2013, MNRAS, 432, 2298
  * Bekki (2015) Bekki K., 2015, MNRAS, 449, 1625
  * Bekki (2017a) Bekki K., 2017a, MNRAS, 467, 1857
  * Bekki (2017b) Bekki K., 2017b, MNRAS, 469, 2933
  * Bekki (2019) Bekki K., 2019, MNRAS, 486, 2570
  * Bekki & Mackey (2009) Bekki K., Mackey A. D., 2009, MNRAS, 394, 124
  * Bekki et al. (2017) Bekki K., Jeřábková T., Kroupa P., 2017, MNRAS, 471, 2242
  * Bellini et al. (2009) Bellini A., Piotto G., Bedin L. R., King I. R., Anderson J., Milone A. P., Momany Y., 2009, A&A, 507, 1393
  * Bellini et al. (2015) Bellini A., et al., 2015, ApJ, 810, L13
  * Bellini et al. (2018) Bellini A., et al., 2018, ApJ, 853, 86
  * Belokurov et al. (2018) Belokurov V., Erkal D., Evans N. W., Koposov S. E., Deason A. J., 2018, MNRAS, 478, 611
  * Bianchini et al. (2013) Bianchini P., Varri A. L., Bertin G., Zocchi A., 2013, ApJ, 772, 67
  * Bianchini et al. (2018) Bianchini P., van der Marel R. P., del Pino A., Watkins L. L., Bellini A., Fardal M. A., Libralato M., Sills A., 2018, MNRAS, 481, 2125
  * Bondi (1952) Bondi H., 1952, MNRAS, 112, 195
  * Bresolin (2019) Bresolin F., 2019, MNRAS, 488, 3826
  * Briley (1997) Briley M. M., 1997, AJ, 114, 1051
  * Calura et al. (2019) Calura F., D’Ercole A., Vesperini E., Vanzella E., Sollima A., 2019, MNRAS, 489, 3269
  * Cameron & Fowler (1971) Cameron A. G. W., Fowler W. A., 1971, ApJ, 164, 111
  * Carretta (2019) Carretta E., 2019, A&A, 624, A24
  * Carretta et al. (2009a) Carretta E., et al., 2009a, A&A, 505, 117
  * Carretta et al. (2009b) Carretta E., Bragaglia A., Gratton R., D’Orazi V., Lucatello S., 2009b, A&A, 508, 695
  * Carretta et al. (2010a) Carretta E., Bragaglia A., D’Orazi V., Lucatello S., Gratton R. G., 2010a, A&A, 519, A71
  * Carretta et al. (2010b) Carretta E., et al., 2010b, A&A, 520, A95
  * Carretta et al. (2011) Carretta E., Lucatello S., Gratton R. G., Bragaglia A., D’Orazi V., 2011, A&A, 533, A69
  * Chantereau et al. (2020) Chantereau W., Biernacki P., Martig M., Bastian N., Salaris M., Teyssier R., 2020, MNRAS, 493, 1306
  * Choi & Yi (2008) Choi E., Yi S. K., 2008, MNRAS, 386, 1332
  * Conroy & Spergel (2011) Conroy C., Spergel D. N., 2011, ApJ, 726, 36
  * Cordoni et al. (2020) Cordoni G., Milone A. P., Mastrobuono-Battisti A., Marino A. F., Lagioia E. P., Tailo M., Baumgardt H., Hilker M., 2020, ApJ, 889, 18
  * Cottrell & Da Costa (1981) Cottrell P. L., Da Costa G. S., 1981, ApJ, 245, L79
  * D’Antona et al. (2012) D’Antona F., D’Ercole A., Carini R., Vesperini E., Ventura P., 2012, MNRAS, 426, 1710
  * D’Antona et al. (2016) D’Antona F., Vesperini E., D’Ercole A., Ventura P., Milone A. P., Marino A. F., Tailo M., 2016, MNRAS, 458, 2122
  * D’Ercole et al. (2008) D’Ercole A., Vesperini E., D’Antona F., McMillan S. L. W., Recchi S., 2008, MNRAS, 391, 825
  * D’Ercole et al. (2012) D’Ercole A., D’Antona F., Carini R., Vesperini E., Ventura P., 2012, MNRAS, 423, 1521
  * D’Orazi & Marino (2010) D’Orazi V., Marino A. F., 2010, ApJ, 716, L166
  * Daddi et al. (2010) Daddi E., et al., 2010, ApJ, 713, 686
  * Dalessandro et al. (2018) Dalessandro E., et al., 2018, ApJ, 859, 15
  * Decressin et al. (2007) Decressin T., Meynet G., Charbonnel C., Prantzos N., Ekström S., 2007, A&A, 464, 1029
  * Decressin et al. (2008) Decressin T., Baumgardt H., Kroupa P., 2008, A&A, 492, 101
  * Denissenkov & Hartwick (2014) Denissenkov P. A., Hartwick F. D. A., 2014, MNRAS, 437, L21
  * Di Criscienzo et al. (2010) Di Criscienzo M., Ventura P., D’Antona F., Milone A., Piotto G., 2010, Monthly Notices of the Royal Astronomical Society, 408, 999
  * Doherty et al. (2014) Doherty C. L., Gil-Pons P., Lau H. H. B., Lattanzio J. C., Siess L., Campbell S. W., 2014, MNRAS, 441, 582
  * Feast et al. (2010) Feast M. W., Abedigamba O. P., Whitelock P. A., 2010, MNRAS, 408, L76
  * Ferraro et al. (2009) Ferraro F. R., et al., 2009, Nature, 462, 483
  * Forbes (2020) Forbes D. A., 2020, MNRAS, 493, 847
  * Fu et al. (2018) Fu X., Bressan A., Marigo P., Girardi L., Montalbán J., Chen Y., Nanni A., 2018, MNRAS, 476, 496
  * Fukui et al. (2017) Fukui Y., Tsuge K., Sano H., Bekki K., Yozin C., Tachihara K., Inoue T., 2017, PASJ, 69, L5
  * Gerber et al. (2020) Gerber J. M., Friel E. D., Vesperini E., 2020, AJ, 159, 50
  * Gratton et al. (2012) Gratton R. G., Carretta E., Bragaglia A., 2012, A&ARv, 20, 50
  * Gratton et al. (2019) Gratton R., Bragaglia A., Carretta E., D’Orazi V., Lucatello S., Sollima A., 2019, A&ARv, 27, 8
  * Harris (1996) Harris W. E., 1996, AJ, 112, 1487
  * Helmi et al. (2018) Helmi A., Babusiaux C., Koppelman H. H., Massari D., Veljanoski J., Brown A. G. A., 2018, Nature, 563, 85
  * Hénault-Brunet et al. (2015) Hénault-Brunet V., Gieles M., Agertz O., Read J. I., 2015, MNRAS, 450, 1164
  * Horta et al. (2020) Horta D., et al., 2020, MNRAS,
  * Hunter (2007) Hunter J. D., 2007, Computing in Science & Engineering, 9, 90
  * Husser et al. (2020) Husser T.-O., et al., 2020, arXiv e-prints, p. arXiv:2001.07725
  * Iben (1975) Iben I. J., 1975, ApJ, 196, 525
  * Iben (1976) Iben I. J., 1976, ApJ, 208, 165
  * Kamann et al. (2018) Kamann S., et al., 2018, MNRAS, 473, 5591
  * Kamann et al. (2020) Kamann S., et al., 2020, MNRAS, 492, 966
  * Karakas (2010) Karakas A. I., 2010, MNRAS, 403, 1413
  * Karakas & Lattanzio (2014) Karakas A. I., Lattanzio J. C., 2014, Publ. Astron. Soc. Australia, 31, e030
  * Kroupa (2019) Kroupa P., 2019, arXiv e-prints, p. arXiv:1910.06971
  * Kruijssen et al. (2020) Kruijssen J. M. D., et al., 2020, MNRAS, 498, 2472
  * Kuzma et al. (2016) Kuzma P. B., Da Costa G. S., Mackey A. D., Roderick T. A., 2016, Monthly Notices of the Royal Astronomical Society, 461, 3639
  * Kuzma et al. (2018) Kuzma P. B., Da Costa G. S., Mackey A. D., 2018, MNRAS, 473, 2881
  * Lagioia et al. (2019) Lagioia E. P., Milone A. P., Marino A. F., Cordoni G., Tailo M., 2019, AJ, 158, 202
  * Lardo et al. (2011) Lardo C., Bellazzini M., Pancino E., Carretta E., Bragaglia A., Dalessandro E., 2011, A&A, 525, A114
  * Larsen et al. (2015) Larsen S. S., Baumgardt H., Bastian N., Brodie J. P., Grundahl F., Strader J., 2015, ApJ, 804, 71
  * Latour et al. (2018) Latour M., Randall S. K., Calamida A., Geier S., Moehler S., 2018, A&A, 618, A15
  * Latour et al. (2019) Latour M., et al., 2019, Astronomy & Astrophysics, 631, A14
  * Layden & Sarajedini (2000) Layden A. C., Sarajedini A., 2000, AJ, 119, 1760
  * Lee & Sneden (2020) Lee J.-W., Sneden C., 2020, arXiv e-prints, p. arXiv:2006.01274
  * Lee et al. (1999) Lee Y. W., Joo J. M., Sohn Y. J., Rey S. C., Lee H. C., Walker A. R., 1999, Nature, 402, 55
  * Magrini et al. (2016) Magrini L., Coccato L., Stanghellini L., Casasola V., Galli D., 2016, A&A, 588, A91
  * Marino et al. (2019) Marino A. F., et al., 2019, Monthly Notices of the Royal Astronomical Society, 487, 3815
  * Marshall et al. (2004) Marshall J. R., van Loon J. T., Matsuura M., Wood P. R., Zijlstra A. A., Whitelock P. A., 2004, MNRAS, 355, 1348
  * Martell (2011) Martell S. L., 2011, Reviews in Modern Astronomy, 23, 173
  * Martell et al. (2016) Martell S. L., et al., 2016, ApJ, 825, 146
  * Mastrobuono-Battisti & Perets (2016) Mastrobuono-Battisti A., Perets H. B., 2016, ApJ, 823, 61
  * McKenzie & Bekki (2018) McKenzie M., Bekki K., 2018, MNRAS, 479, 3126
  * McKinney (2010) McKinney 2010, in Stéfan van der Walt Jarrod Millman eds, Proceedings of the 9th Python in Science Conference. pp 56 – 61, doi:10.25080/Majora-92bf1922-00a
  * Milone (2019) Milone A. P., 2019, arXiv e-prints, p. arXiv:1908.11703
  * Milone et al. (2012) Milone A. P., et al., 2012, ApJ, 744, 58
  * Milone et al. (2015) Milone A. P., et al., 2015, ApJ, 808, 51
  * Milone et al. (2017) Milone A. P., et al., 2017, MNRAS, 464, 3636
  * Milone et al. (2018) Milone A. P., Marino A. F., Mastrobuono-Battisti A., Lagioia E. P., 2018, MNRAS, 479, 5005
  * Milone et al. (2020) Milone A. P., et al., 2020, MNRAS, 491, 515
  * Myeong et al. (2018) Myeong G. C., Evans N. W., Belokurov V., Sand ers J. L., Koposov S. E., 2018, ApJ, 863, L28
  * Myeong et al. (2019) Myeong G. C., Vasiliev E., Iorio G., Evans N. W., Belokurov V., 2019, MNRAS, 488, 1235
  * Mészáros et al. (2019) Mészáros S., et al., 2019, Monthly Notices of the Royal Astronomical Society, 492, 1641
  * Nardiello et al. (2018) Nardiello D., et al., 2018, MNRAS, 477, 2004
  * Navarro et al. (1996) Navarro J. F., Frenk C. S., White S. D. M., 1996, ApJ, 462, 563
  * Neto et al. (2007) Neto A. F., et al., 2007, MNRAS, 381, 1450
  * Niederhofer et al. (2017) Niederhofer F., et al., 2017, MNRAS, 465, 4159
  * Norris & Freeman (1979) Norris J., Freeman K. C., 1979, ApJ, 230, L179
  * Osman et al. (2020) Osman O., Bekki K., Cortese L., 2020, arXiv e-prints, p. arXiv:2009.08078
  * Pancino et al. (2017) Pancino E., et al., 2017, Astronomy & Astrophysics, 601, A112
  * Pérez & Granger (2007) Pérez F., Granger B. E., 2007, Computing in Science and Engineering, 9, 21
  * Pflamm-Altenburg & Kroupa (2009) Pflamm-Altenburg J., Kroupa P., 2009, MNRAS, 397, 488
  * Piotto et al. (2015) Piotto G., et al., 2015, The Astronomical Journal, 149, 91
  * Prantzos & Charbonnel (2006) Prantzos N., Charbonnel C., 2006, A&A, 458, 135
  * Renzini (2008) Renzini A., 2008, MNRAS, 391, 354
  * Renzini et al. (2015) Renzini A., et al., 2015, MNRAS, 454, 4197
  * Rosen & Bregman (1995) Rosen A., Bregman J. N., 1995, ApJ, 440, 634
  * Rossi et al. (2016) Rossi L. J., Bekki K., Hurley J. R., 2016, MNRAS, 462, 2861
  * Sbordone et al. (2011) Sbordone L., Salaris M., Weiss A., Cassisi S., 2011, A&A, 534, A9
  * Schiavon et al. (2016) Schiavon R. P., et al., 2016, Monthly Notices of the Royal Astronomical Society, 465, 501
  * Searle (1971) Searle L., 1971, ApJ, 168, 327
  * Siess (2010) Siess L., 2010, A&A, 512, A10
  * Smith & Norris (1982) Smith G. H., Norris J., 1982, ApJ, 254, 149
  * Sollima et al. (2019) Sollima A., Baumgardt H., Hilker M., 2019, MNRAS, 485, 1460
  * Sutherland & Dopita (1993) Sutherland R. S., Dopita M. A., 1993, ApJS, 88, 253
  * Szécsi & Wünsch (2019) Szécsi D., Wünsch R., 2019, ApJ, 871, 20
  * Tsuge et al. (2019) Tsuge K., et al., 2019, ApJ, 871, 44
  * Van Der Walt et al. (2011) Van Der Walt S., Colbert S. C., Varoquaux G., 2011, Computing in Science & Engineering, 13, 22
  * Ventura & D’Antona (2008a) Ventura P., D’Antona F., 2008a, MNRAS, 385, 2034
  * Ventura & D’Antona (2008b) Ventura P., D’Antona F., 2008b, A&A, 479, 805
  * Ventura & D’Antona (2011) Ventura P., D’Antona F., 2011, MNRAS, 410, 2760
  * Vesperini et al. (2010) Vesperini E., McMillan S. L. W., D’Antona F., D’Ercole A., 2010, ApJ, 718, L112
  * Vesperini et al. (2018) Vesperini E., Hong J., Webb J. J., D’Antona F., D’Ercole A., 2018, MNRAS, 476, 2731
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Vorobyov et al. (2004) Vorobyov E. I., Klein U., Shchekinov Y. A., Ott J., 2004, A&A, 413, 939
  * Walter et al. (2008) Walter F., Brinks E., de Blok W. J. G., Bigiel F., Kennicutt Robert C. J., Thornley M. D., Leroy A., 2008, AJ, 136, 2563
  * Wang et al. (2020) Wang L., Kroupa P., Takahashi K., Jerabkova T., 2020, MNRAS, 491, 440
  * Warren et al. (2011) Warren S. R., et al., 2011, The Astrophysical Journal, 738, 10
  * de Mink et al. (2009) de Mink S. E., Pols O. R., Langer N., Izzard R. G., 2009, A&A, 507, L1
  * van den Bergh (1996) van den Bergh S., 1996, ApJ, 471, L31

##  Appendix A Spatial distributions for the remaining components within the
simulation

The companion figures of Fig 3 for the gas and disc star components.

Figure 15: The spatial distribution of gas particles during the simulation.
This includes both ISM and AGB ejecta. The figure shares the same properties
as Fig. 3.  Figure 16: The spatial distribution of disc stars during the
simulation. The figure shares the same properties as Fig. 3.

##  Appendix B Simulation Models

We present the results of all 370 Myr simulations which were used to construct
the scaling relation in Fig 12. We do not include our larger sample of 140 Myr
simulations which were used to inform our current parameter choices.

Table 4: Parameters for our Galactic GC models which ran for 370 Myr. $\rm{R}_{start}$ denotes the distance of the GCs starting location with respect to the host galaxy’s centre of mass. Model ID |  $M_{\rm{dm}}$ ($\rm{M}_{\odot}$) |  $\rm{R}_{galaxy}$ (kpc) |  $M_{gc}$ ($\rm{M}_{\odot}$) |  $\rm{R}_{gc}$ (pc) |  $\rm{R}_{HI\ hole}$ (kpc) | $f_{\rm rot}$ |  $\theta_{\rm gc}$ (∘) |  $\rm{R}_{start}$ (pc)
---|---|---|---|---|---|---|---|---
M1 | 5$\times 10^{9}$ | 2.4 | 1$\times 10^{6}$ | 0.02 | 0.2 | 0.2 | 0 | 85
M2 | 5$\times 10^{9}$ | 1.7 | 1$\times 10^{6}$ | 0.02 | 0.1 | 0.2 | 30 | 36
M3 | 5$\times 10^{9}$ | 1.7 | 1$\times 10^{6}$ | 0.02 | 0.1 | 0.2 | 30 | 187
M4 | 5$\times 10^{9}$ | 1.7 | 1$\times 10^{6}$ | 0.02 | 0.3 | 0.2 | 30 | 36
M5 | 5$\times 10^{9}$ | 1.7 | 1$\times 10^{6}$ | 0.02 | 0.3 | 0.2 | 30 | 36
M6 | 5$\times 10^{9}$ | 1.7 | 1$\times 10^{6}$ | 0.02 | 0.2 | 0.2 | 30 | 179
M7 | 5$\times 10^{9}$ | 1.7 | 3$\times 10^{5}$ | 0.02 | 0.3 | 0.2 | 30 | 36
M8 | 5$\times 10^{9}$ | 1.7 | 3$\times 10^{5}$ | 0.02 | 0.2 | 0.2 | 30 | 36
M9 | 5$\times 10^{9}$ | 1.7 | 1$\times 10^{6}$ | 0.02 | 0.2 | 0.2 | 30 | 38
M10 | 5$\times 10^{9}$ | 1.7 | 1$\times 10^{6}$ | 0.02 | 0.2 | 0.2 | 0 | 130
M11 | 5$\times 10^{9}$ | 1.7 | 1$\times 10^{6}$ | 0.02 | 0.2 | 0.2 | 0 | 370
M12 | 5$\times 10^{9}$ | 1.7 | 3$\times 10^{5}$ | 0.02 | 0.2 | 0.2 | 30 | 36
M13 | 5$\times 10^{9}$ | 1.7 | 1$\times 10^{6}$ | 0.02 | 0.2 | 0.2 | 0 | 83
M14 | 5$\times 10^{9}$ | 1.7 | 1$\times 10^{6}$ | 0.02 | 0.2 | 0.2 | 0 | 466
M15 | 5$\times 10^{9}$ | 2.4 | 1$\times 10^{6}$ | 0.02 | 0.1 | 0.2 | 0 | 85
M16 | 5$\times 10^{9}$ | 2.4 | 1$\times 10^{6}$ | 0.02 | 0.2 | 0.2 | 0 | 450
M17 | 2$\times 10^{10}$ | 3.4 | 1$\times 10^{6}$ | 0.02 | 0.2 | 0.2 | 0 | 140
M18 | 5$\times 10^{9}$ | 1.7 | 1$\times 10^{6}$ | 0.02 | 0.2 | 0.2 | 0 | 152
M19 | 2$\times 10^{10}$ | 4.9 | 1$\times 10^{6}$ | 0.02 | 0.2 | 0.2 | 0 | 152
M20 | 5$\times 10^{9}$ | 2.4 | 1$\times 10^{6}$ | 0.02 | 0.2 | 0.2 | 0 | 18
M21 | 5$\times 10^{9}$ | 2.4 | 1$\times 10^{5}$ | 0.02 | 0.2 | 0.2 | 0 | 173
M22 | 5$\times 10^{9}$ | 2.4 | 5$\times 10^{5}$ | 0.02 | 0.2 | 0.2 | 0 | 159
M23 | 5$\times 10^{9}$ | 2.4 | 1$\times 10^{6}$ | 0.02 | 0.2 | 0.2 | 0 | 580
M24 | 5$\times 10^{9}$ | 2.4 | 1$\times 10^{6}$ | 0.02 | 0.2 | 0.03 | 0 | 159
M25 | 5$\times 10^{9}$ | 2.4 | 5$\times 10^{5}$ | 0.02 | 0.2 | 0.03 | 0 | 159
M26 | 5$\times 10^{9}$ | 2.4 | 1$\times 10^{6}$ | 0.02 | 0.2 | 0.03 | 0 | 991
M27 | 5$\times 10^{9}$ | 2.4 | 1$\times 10^{6}$ | 0.02 | 0.2 | 0.2 | 60 | 105
M28 | 5$\times 10^{9}$ | 1.7 | 1$\times 10^{6}$ | 0.02 | 0.1 | 0.2 | 30 | 141
M29 | 5$\times 10^{9}$ | 2.4 | 1$\times 10^{6}$ | 0.02 | 0.2 | 0.2 | 0 | 78
M30 | 5$\times 10^{9}$ | 2.4 | 1$\times 10^{6}$ | 0.02 | 0.2 | 0.2 | 0 | 81
M31 | 5$\times 10^{9}$ | 2.4 | 1$\times 10^{6}$ | 0.02 | 0.2 | 0.2 | 0 | 81
M32 | 5$\times 10^{9}$ | 2.4 | 1$\times 10^{6}$ | 0.02 | 0.2 | 0.2 | 0 | 79
M33 | 5$\times 10^{9}$ | 2.4 | 3$\times 10^{5}$ | 0.02 | 0.2 | 0.2 | 0 | 108
M34 | 5$\times 10^{9}$ | 3.4 | 1$\times 10^{6}$ | 0.02 | 0.2 | 0.2 | 0 | 86
M35 | 5$\times 10^{9}$ | 2.4 | 3$\times 10^{6}$ | 0.02 | 0.2 | 0.2 | 0 | 385
Table 5: Parameters for our AGB models. See the Galactic GC model table for a description of the headers. Model ID |  $M_{\rm{dm}}$ ($\rm{M}_{\odot}$) |  $\rm{R}_{galaxy}$ (kpc) |  $M_{gc}$ ($\rm{M}_{\odot}$) |  $\rm{R}_{gc}$ (pc) |  $\rm{R}_{HI\ hole}$ (kpc) | $f_{\rm rot}$ |  $\theta_{\rm gc}$ (∘) |  $\rm{R}_{start}$ (pc)
---|---|---|---|---|---|---|---|---
M36 | 5$\times 10^{9}$ | 1.7 | 1$\times 10^{6}$ | 0.02 | 0.5 | 0.2 | 30 | 36
M37 | 5$\times 10^{9}$ | 2.4 | 1$\times 10^{6}$ | 0.02 | 0.5 | 0.03 | 0 | 241
M38 | 5$\times 10^{9}$ | 2.4 | 5$\times 10^{6}$ | 0.02 | 0.5 | 0.03 | 0 | 223
M39 | 5$\times 10^{9}$ | 2.4 | 5$\times 10^{6}$ | 0.02 | 0.5 | 0.2 | 0 | 223
M40 | 5$\times 10^{9}$ | 2.4 | 1$\times 10^{6}$ | 0.02 | 0.5 | 0.2 | 60 | 402
M41 | 5$\times 10^{9}$ | 2.4 | 1$\times 10^{6}$ | 0.02 | 1.0 | 0.2 | 0 | 85
M42 | 5$\times 10^{9}$ | 2.4 | 1$\times 10^{6}$ | 0.02 | 2.0 | 0.2 | 0 | 79
M43 | 5$\times 10^{9}$ | 2.4 | 3$\times 10^{6}$ | 0.02 | 1.0 | 0.2 | 0 | 383
Table 6: Parameters for our LSB galaxy models. See the Galactic GC model table for a description of the headers. Model ID |  $M_{\rm{dm}}$ ($\rm{M}_{\odot}$) |  $\rm{R}_{galaxy}$ (kpc) |  $M_{gc}$ ($\rm{M}_{\odot}$) |  $\rm{R}_{gc}$ (pc) |  $\rm{R}_{HI\ hole}$ (kpc) | $f_{\rm rot}$ |  $\theta_{\rm gc}$ (∘) |  $\rm{R}_{start}$ (pc)
---|---|---|---|---|---|---|---|---
M44 | 5$\times 10^{9}$ | 2.4 | 3$\times 10^{5}$ | 0.02 | 0.2 | 0.2 | 0 | 1907
M45 | 5$\times 10^{9}$ | 4.9 | 1$\times 10^{6}$ | 0.02 | 0.2 | 0.2 | 0 | 196
M46 | 5$\times 10^{9}$ | 2.4 | 1$\times 10^{5}$ | 0.02 | 0.2 | 0.2 | 0 | 1941
M47 | 5$\times 10^{9}$ | 2.4 | 1$\times 10^{6}$ | 0.02 | 0.2 | 0.2 | 0 | 1951