Massive star clusters in the semi-analytical galaxy formation model L-Galaxies 2020

Nils Hoyer; Silvia Bonoli; Nate Bastian; Diego Herrero-Carrión; Nadine Neumayer; David Izquierdo-Villalba; Daniele Spinoso; Robert M. Yates; Markos Polkas; Maria Celeste Artale

doi:10.1051/0004-6361/202554325

Home

All issues

Volume 709 (May 2026)

A&A, 709 (2026) A109

Full HTML

Open Access

Issue		A&A Volume 709, May 2026


Article Number		A109
Number of page(s)		24
Section		Galactic structure, stellar clusters and populations
DOI		https://doi.org/10.1051/0004-6361/202554325
Published online		07 May 2026

A&A, 709, A109 (2026)

Massive star clusters in the semi-analytical galaxy formation model L-Galaxies 2020

Nils Hoyer¹^,2^,3^,4^★, Silvia Bonoli²^,5, Nate Bastian²^,5, Diego Herrero-Carrión²^,6, Nadine Neumayer³, David Izquierdo-Villalba⁷, Daniele Spinoso⁸, Robert M. Yates⁹, Markos Polkas¹⁰^,2^,5 and Maria Celeste Artale¹¹

¹ LIRA, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, Université Paris Cité, CY Cergy Paris Université, 5 place Jules Janssen, 92195 Meudon, France
² Donostia International Physics Center, Paseo Manuel de Lardizabal 4, 20118 Donosita-San Sebastián, Spain
³ Max-Planck-Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany
⁴ Universität Heidelberg, Seminarstrasse 2, 69117 Heidelberg, Germany
⁵ Ikerbasque, Basque Foundation for Science, 48013 Bilbao, Spain
⁶ University of the Basque Country UPV/EHU, Department of Theoretical Physics, Bilbao 48080, Spain
⁷ Institute of Space Sciences (ICE, CSIC), Campus UAB, Carrer de Magrans, 08193 Barcelona, Spain
⁸ Dipartimento di Fisica, Universitá degli Studi di Milano Bicocca, piazza della Scienza 3, 20126 Milano, Italy
⁹ Centre for Astrophysics Research, University of Hertfordshire, Hatfield, AL10 9AB, UK
¹⁰ School of Physics and Astronomy & Institute for Gravitational Wave Astronomy, University of Birmingham, Birmingham, B15 2TT, UK
¹¹ Departamento de Ciencias Físicas, Facultad de Ciencias Exactas, Universidad Andres Bello, Fernández Concha 700, Las Condes, Santiago, Chile

^★ Corresponding author: This email address is being protected from spambots. You need JavaScript enabled to view it.

Received: 28 February 2025
Accepted: 23 March 2026

Abstract

It has been established that a direct link exists between the formation history of star cluster populations and their host galaxies. However, our lack of understanding of the assembly of star cluster populations impedes our ability to use them as tracers of galaxy evolution. In this work we introduce a new variation of the L-Galaxies 2020 semi-analytic galaxy formation model that includes the formation of star clusters above 10⁴ M_⊙ and probes different physical assumptions that affect their evolution over cosmic time. We used properties of different galaxy components and localised star formation to determine the bound fraction of star formation in disks. After randomly sampling masses from an environmentally dependent star cluster initial mass function, we assigned to each object a half-mass radius, metallicity, and distance from the galaxy centre. We considered up to 2000 individual star clusters per galaxy and evolved their properties over time while taking into account stellar evolution, two-body relaxation, tidal shocks, dynamical friction, and re-positioning during galaxy mergers. Our simulation successfully reproduces several observational quantities, such as the empirical relationship between the absolute V-band magnitude of the brightest young star clusters and the host galaxy star formation rate, the mass function of young star clusters, and the mean metallicities of the star cluster distributions versus galaxy masses. The simulation reveals great complexity in the z = 0 star cluster population resulting from differential destruction channels and origins, including in situ populations in the disk, a major merger-induced heated component in the halo, and accreted star clusters. Model variations point out the importance of the shape of the star cluster initial mass function, the initial distribution of half-mass radii, and the relationship between the sound speed of cold gas and the star formation rate. Our new model provides new avenues to trace individual star clusters and test cluster-related physics within a cosmological set-up in a computationally efficient manner.

Key words: galaxies: star clusters: general

© The Authors 2026

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

This article is published in open access under the Subscribe to Open model. This email address is being protected from spambots. You need JavaScript enabled to view it. to support open access publication.

1 Introduction

A natural consequence of star formation in cold gas-dense and rotationally unstable regions is the formation of bound stellar structures, ranging from few-body systems to massive star clusters whose mass rivals the baryonic components of entire (dwarf) galaxies. The properties and survival times of these star clusters heavily depend on their changing environment, which makes them excellent tracers of galaxy assembly. To date, a large body of work has used the present-day star cluster population, such as globular clusters, to constrain the evolutionary history of their current host galaxy, including the Milky Way, M 31, and other nearby systems (e.g. Huchra et al. 1991; Barmby et al. 2000; Perrett et al. 2002; West et al. 2004; Brodie & Strader 2006; Forbes & Bridges 2010; Leaman et al. 2013; Huxor et al. 2014; Cantiello et al. 2015; Veljanoski et al. 2015; Myeong et al. 2018; Callingham et al. 2022; Hammer et al. 2023, 2024; Ines Ennis et al. 2024; Usher et al. 2024). However, despite the clear link between galaxy and star cluster formation that is expressed via, for example, the scaling relation between the number of star clusters and host galaxy mass (e.g. West et al. 1995; Blakeslee 1999; Peng et al. 2008; Spitler & Forbes 2009; Hudson et al. 2014; El-Badry et al. 2019; Bastian et al. 2020; Burkert & Forbes 2020; Zaritsky 2022; Le & Cooper 2025), precise details about the formation environments and initial properties of z = 0 star clusters remain elusive (e.g. Forbes et al. 2018; Valenzuela et al. 2025).

Constraints on the ages of globular clusters in galaxies, including the Milky Way, suggest that they form at redshifts z ≳ 1 (e.g. Carretta et al. 2000; Krauss & Chaboyer 2003; Lee 2003; Kaviraj et al. 2005; Strader et al. 2005; Correnti et al. 2016). Studying the natal environment of clusters at these distances is challenging; however, observations with the Very Large Telescope, the Hubble Space Telescope, and the James Webb Space Telescope, combined with strong gravitational lensing, make such observations possible (e.g. Vanzella et al. 2017, 2019; Kikuchihara et al. 2020; Mowla et al. 2022; Forbes & Romanowsky 2023; Vanzella et al. 2023; Claeyssens et al. 2023; Mowla et al. 2024). Most recently, Adamo et al. (2024) discussed the properties of young and massive star clusters in the Cosmic Gems arc at redshift z ≈ 10 (Salmon et al. 2018; Bradley et al. 2024). Such clusters are compact (half-light radius ≤2 pc), massive (stellar mass ≳10⁶ M_⊙), and could potentially contribute to the z = 0 star cluster population. Irrespective of whether they survive until z = 0, they constitute a significant baryonic component of high-z galaxies (≳30 % for the Cosmic Gems Arc; Adamo et al. 2024) and most likely influence the host galaxy’s evolution.

Simulations are required to understand the properties of the full star cluster distribution, as observations can only trace the brightest and/or most massive and most unobscured star clusters. Detailed simulations of star cluster formation and evolution have now become feasible in high-resolution hydrodynamical simulations of galaxy formation, including from dwarf galaxies to Milky Way analogues. These simulations either include star clusters in full cosmological simulations, which sometimes include adaptive mesh refinement or zoom-in techniques (e.g. Li et al. 2017, 2018; Li & Gnedin 2019; Brown & Gnedin 2022; Reina-Campos et al. 2022a; Garcia et al. 2023; Calura et al. 2025), or use high spatial and temporal resolutions but evolve for less than one giga year (e.g. Lahén et al. 2020; Hislop et al. 2022; Lahén et al. 2023; Elmegreen & Lahén 2024; Lahén et al. 2024, 2025; Reina-Campos et al. 2025). Another approach is to add star cluster-related physics to existing simulations. For example, this approach was realised in the E-MOSAICS project (Pfeffer et al. 2018; Kruijssen et al. 2019a), where prescriptions for star cluster formation were added in post-processing to the EAGLE simulation (Crain et al. 2015; Schaye et al. 2015). Focussing on E-MOSAICS, their simulation can reproduce a number of observables related to the masses and metallicity distributions of the globular clusters as well as relations with properties of their host galaxy and dark matter halo (see details in Hughes et al. 2019; Kruijssen et al. 2019a,b; Reina-Campos et al. 2019; Pfeffer et al. 2019a,b; Bastian et al. 2020; Hughes et al. 2020; Kruijssen et al. 2020; Reina-Campos et al. 2020; Keller et al. 2020; Horta et al. 2021; Trujillo-Gomez et al. 2021; Dolfi et al. 2022; Hughes et al. 2022; Reina-Campos et al. 2022b, 2023b; Pfeffer et al. 2024; Newton et al. 2025; Pfeffer et al. 2025).

All of the above mentioned simulations follow in detail the evolution of stars within individual star clusters or of individual baryonic particles, which contain star clusters but are expensive to run. In contrast, semi-analytical models sacrifice resolution in order to lower the computational expense and gain the ability to explore a wide parameter space of various astrophysical mechanisms (e.g. White & Rees 1978; White & Frenk 1991; Baugh 2006; Somerville et al. 2008; Somerville & Davé 2015, but see also Valenzuela et al. 2021 for an empirical model). This approach has proven to be successful in reproducing observational quantities related to, for example, the co-evolution of galaxies and massive black holes (e.g. Kauffmann & Haehnelt 2000; Croton et al. 2006; Croton 2006; De Lucia & Blaizot 2007; Monaco et al. 2007; Bonoli et al. 2009, 2014; Izquierdo-Villalba et al. 2020; Gabrielpillai et al. 2022; del P. Lagos et al. 2024), which allows the values of free parameters of the assumed physical models to be constrained.

For star clusters in a semi-analytical galaxy formation framework specifically, most recently De Lucia et al. (2024) added prescriptions for their formation to the GAEA model (De Lucia et al. 2014; Hirschmann et al. 2016). The authors were able to reproduce the empirical relationship between the total mass in globular clusters and the parent halo mass (see also e.g. Kravtsov & Gnedin 2005; Prieto & Gnedin 2008; Muratov & Gnedin 2010; Li et al. 2017) but utilised simple prescriptions related to mass loss and dynamical friction and relied on global star cluster population statistics.

In this first paper of a planned series, we introduce an implementation of the formation and evolution of star clusters into a public version of the semi-analytical galaxy formation model ‘L-Galaxies’ (Henriques et al. 2020; Yates et al. 2021). Our work differs from the above mentioned work by De Lucia et al. (2024) in that we track the evolution of individual clusters and use different sets of astrophysical prescriptions that make use of the radially resolved gas and stellar disks in the 2020 version of the code. This effort enabled us to study individual star clusters across different galaxy types, masses, environments, and redshifts, and it offers new avenues to study the formation of nuclear star clusters and the co-evolution of black holes with star clusters in future work.

We begin in Section 2 by detailing the governing equations of the model, starting with galaxy components to evaluate the formation efficiency of star clusters, the initial properties of the star clusters, and eventually the evolution of star clusters within the evolution of their host galaxies. We then evaluate the results of our model in Section 3, focussing on young massive star clusters in disk-dominated galaxies and the metallicity distributions of in situ and accreted star clusters for different galaxy morphologies. We also discuss the caveats of our approach. We conclude in Section 4 and present an outlook for future papers in this series. Appendix A presents variations of key model parameters.

2 Model description

Below we outline the basic principles behind our model that we summarise in Figure 1. Going forward, when mentioning the term L-Galaxies we specifically refer to a modified version of the 2020 model introduced by Yates et al. (2021). This version of the code improves over the default model (Henriques et al. 2020) by modifying the prescriptions for metal injection from stellar winds and supernovae into the circum-galactic medium.

2.1 The L-Galaxies semi-analytical galaxy formation model

The L-Galaxies model combines merger trees from dark matter-only N-body simulations with a set of partial differential equations for the evolution of baryonic components. It has been developed to primarily run on the Millennium (Springel et al. 2005) and Millennium-II (Boylan-Kolchin et al. 2009) simulations with box sizes and dark matter particle masses of 480.3 h⁻¹Mpc and 9.61 × 10⁸ h⁻¹M_⊙, as well as 96.1 h⁻¹Mpc and 7.69 × 10⁶ h⁻¹M_⊙, respectively. Dark matter (sub-)halos are identified using a ‘friends-of-friends’ algorithm (Davis et al. 1985) and the ‘subfind’ algorithm (Springel et al. 2001; Dolag et al. 2009) and are used as input to L-Galaxies. As discussed in the works related to the last few major releases (Guo et al. 2011; Henriques et al. 2015, 2020), the model can reproduce many observables of baryonic components, such as the redshift-dependent galaxy mass function, passive galaxy fraction, and the cosmic density of the star formation rate (SFR).

Many extensions to L-Galaxies have been developed, including those that focus on the gas (Vijayan et al. 2019; Ayromlou et al. 2021; Yates et al. 2021; Parente et al. 2023; Zhong et al. 2023a,b; Parente et al. 2024; Yates et al. 2024), stars (Bluck et al. 2016; Wang et al. 2018; Irodotou et al. 2019; Izquierdo-Villalba et al. 2019; Murphy et al. 2022; Wang & Peng 2024), massive black holes (Bonoli et al. 2009, 2014; Izquierdo-Villalba et al. 2020, 2022; Spinoso et al. 2023; Polkas et al. 2024), and other components (Barrera et al. 2023; Vani et al. 2025). This shows that L-Galaxies is a versatile utility to explore the assembly of galaxies and has the potential to investigate the formation of star clusters as well.

One of the key features of the 2020 version of L-Galaxies that was adapted from Fu et al. (2013) is the introduction of concentric annuli. By default, the model features twelve such annuli that are logarithmically spaced and act as the resolution limit for the cold gas and stars within a galaxy’s disk and stars within a galaxy’s bulge. The annuli’s outer radii have values of $w_{j} / [h^{- 1} kpc] = 0.01 \times 2^{j} with j \in [1, 12],$ $Mathematical equation: ${w_j}/\left[ {{h^{ - 1}}{\rm{kpc}}} \right] = 0.01 \times {2^j}{\rm{with}}j \in \left[ {1,12} \right],$$ (1)

resulting in w₁ ≈ 29.7 pc and w₁₂ ≈ 60.8 kpc.

One of the affected properties of the separation into annuli is the star formation prescription. L-Galaxies assumes that the molecular gas in each annulus collapses on a dynamical timescale, τ_dyn, and is transformed to stars with an efficiency, $\in_{H_{2}}$ $Mathematical equation: ${ \in _{{{\rm{H}}_2}}}$$ (e.g. Fu et al. 2012). Thus, in terms of surface mass density, the SFR for ring j is $Σ_{SFR, j} = ε_{H_{2}} τ_{dyn}^{- 1} Σ_{H_{2}, j} .$ $Mathematical equation: ${\Sigma _{{\rm{SFR}},j}} = {\varepsilon _{{{\rm{H}}_2}}}\tau _{{\rm{dyn}}}^{ - 1}{\Sigma _{{{\rm{H}}_2},j}}.$$ (2)

For the dynamical time, the model assumes $τ_{dyn} = R_{g} / v_{max},$ $Mathematical equation: ${\tau _{{\rm{dyn}}}} = {R_{\rm{g}}}/{v_{{\rm{max}}}},$$ (3)

relating the disk scale length of the cold gas, R_g, and the maximum value of the rotation velocity of the dark matter halo, v_max. The amount of molecular gas itself is derived from the available cold gas in each ring and, in turn, is related to the gas’ metallicity and clumping of gas clouds. Extensive discussions and recipes for computing the molecular mass (per annuli) are presented in Krumholz et al. (2009); Fu et al. (2010); McKee & Krumholz (2010); Fu et al. (2013); Henriques et al. (2020).

In this work we adopted the 2014 cosmology of Planck (Planck Collaboration XVI 2014) with Ω_Λ,0 = 0.685, Ω_M,0 = 0.315 (with Ω_B,0 = 0.0487), σ₈ = 0.826, n_s = 0.96, and h = 0.673. The L-Galaxies itself utilises dark matter-only simulations that are re-scaled (Angulo & White 2010; Angulo & Hilbert 2015) to these values.

Fig. 1

Relevant prescriptions for the assembly and evolution of star clusters implemented into a modified version (Yates et al. 2021) of L-Galaxies 2020 (Henriques et al. 2020).

2.2 Gravitational potential

The computation of the fraction of star formation that is bound in star clusters (Section 2.3.3) requires knowledge of the underlying gravitational potential. We consider the contribution of all available galaxy components that L-Galaxies tracks over time, namely a central massive black hole, dark matter, a bulge, gaseous and stellar disks, and gaseous and stellar halos.

Massive black hole. In case a massive black hole occupies a galaxy centre, we introduced a gravitational potential of a point mass, i.e. $Φ_{BH} (w) = - \frac{G M_{BH}}{w},$ $Mathematical equation: ${\Phi _{{\rm{BH}}}}\left( w \right) = - {{G{M_{{\rm{BH}}}}} \over w},$$ (4)

where G is the gravitational constant, M_BH is the black hole mass, and w is the galactocentric distance. We assume a constant alignment between the dark matter and galaxy centres and that the massive black hole does not ‘wander’ within the galaxy (see Izquierdo-Villalba et al. 2020; Untzaga et al. 2024, for a discussion on wandering black holes in L-Galaxies).

Dark matter halo. We selected the classical Navarro-Frenk-White (NFW) profile (Navarro et al. 1996) to describe the distribution of dark matter. The gravitational potential is given by $Φ_{DM} (w) = - \frac{G M_{vir}}{w} \frac{\ln (1 + c_{vir} w / R_{vir})}{\ln (1 + c_{vir}) - c_{vir} / (1 + c_{vir})},$ $Mathematical equation: ${\Phi _{{\rm{DM}}}}(w) = - {{G{M_{{\rm{vir}}}}} \over w}{{\ln (1 + {c_{{\rm{vir}}}}w/{R_{{\rm{vir}}}})} \over {\ln (1 + {c_{{\rm{vir}}}}) - {c_{{\rm{vir}}}}/(1 + {c_{{\rm{vir}}}})}},$$ (5)

where we introduced the virial mass M_vir and radius R_vir, respectively. For the concentration parameter c_vir, which relates the profile’s scale-radius to R_vir, we assumed a relation from Dutton & Macció (2014) that connects it to M_vir via $\log_{10} c_{vir} = α (z) + β (z) \log_{10} (M_{vir} / 10^{12} h^{- 1} M_{⊙}) .$ $Mathematical equation: ${\log _{10}}{c_{{\rm{vir}}}} = \alpha (z) + \beta (z){\log _{10}}({M_{{\rm{vir}}}}/{10^{12}}{h^{ - 1}}{{\rm{M}}_ \odot }).$$ (6)

The redshift-dependent coefficients α(z) and β(z) take values according to Table 3 of Dutton & Macció (2014) with α(0) = 1.025 and β(0) = −0.097.

Galactic bulge. We modelled a galaxy’s bulge with a Jaffe (1983) profile of the form $Φ_{B} (w) = - \frac{G M_{B}}{w_{B}} \ln (1 + w_{B} / w) .$ $Mathematical equation: ${\Phi _{\rm{B}}}(w) = - {{G{M_{\rm{B}}}} \over {{w_{\rm{B}}}}}\ln (1 + {w_{\rm{B}}}/w).$$ (7)

Its scale length, w_B, encloses half of the bulge’s mass, M_B.

Gaseous and stellar disks. We followed the default assumption in L-Galaxies that both the gaseous and stellar disks are well described by two-dimensional exponential density profiles. The mid-plane gravitational potential of the two disks is expressed with modified Bessel functions of the first and second kind, I_ν and K_ν, respectively (e.g. Watson 1944; Kuijken & Gilmore 1989). The gravitational potential for disk i reads $Φ_{D,i} (w) = - π G Σ_{D,i} w [I_{0} (y_{i}) K_{1} (y_{i}) - I_{1} (y_{i}) K_{0} (y_{i})],$ $Mathematical equation: ${\Phi _{{\rm{D,i}}}}(w) = - \pi G{\Sigma _{{\rm{D,i}}}}w[{I_0}({y_i}){K_1}({y_i}) - {I_1}({y_i}){K_0}({y_i})],$$ (8)

where y_i = w/(2 × w_D,i), w_D,i is a characteristic scale length of disk i (gaseous or stellar), and Σ_D,i is the central surface mass density (Freeman 1970; Binney & Tremaine 2008).

Hot gas and stellar halos. We assumed that a galaxy’s halo contains both hot gas that is unable to form stars and a stellar medium that originates entirely from stripped satellite galaxies. For simplicity, we assumed that an isothermal profile can describe both the gaseous and stellar halo well, i.e. $Φ_{H,i} (w) \propto - \frac{G M_{H,i}}{R_{H,i}} \ln (w) .$ $Mathematical equation: ${\Phi _{{\rm{H,i}}}}(w) \propto - {{G{M_{{\rm{H,i}}}}} \over {{R_{{\rm{H,i}}}}}}\ln (w).$$ (9)

Furthermore, we assumed that all mass is enclosed within the virial radius such that R_H,i = R_vir¹.

2.3 Star cluster formation

The total mass of newly formed stars that are bound in star clusters equals the SFR multiplied by the simulation time step and the cluster formation efficiency (Bastian 2008). In principle, the latter must be computed considering local conditions of the gas phase, mergers from substructures, and accretion of gas during cluster formation (see e.g. Karam & Sills 2022, 2023). This step should ideally be computed at sub-parsec resolution (Renaud 2020) with a star-by-star treatment of feedback (Calura et al. 2025) but this is not feasible with our approach.

Here we follow the prescription outlined by Kruijssen (2012) to estimate the cluster formation efficiency. The model relies on three quantities: (1) the cold gas surface mass density Σ_g, which we compute as an annuli’s molecular mass divided by its area, (2) the epicyclic frequency κ, and (3) the Toomre disk stability parameter Q. We compute these values for the logarithmic mean galactocentric distance of each annuli j and assume that it is applicable to a wide range of environments and redshifts². Furthermore, we assumed that all star clusters form in galaxy disks, thus neglecting cluster formation during galaxy mergers or outside of dark matter halos at high-z (Lake et al. 2021, 2023).

2.3.1 Epicyclic frequency

Based on the combined gravitational potential of all previously introduced galaxy components, ***(11)*** $Φ_{tot} (w) = Σ_{c} Φ_{c} (w)$ $Mathematical equation: ${\Phi _{{\rm{tot}}}}\left( w \right) = {\Sigma _c}{\Phi _c}\left( w \right)$$ , we could easily determine the epicyclic frequency for circular orbits. For annulus j and log-mean distance ⟨wj⟩, $κ_{j} = \sqrt{\frac{3}{〈 w_{j} 〉} {\frac{\partial Φ_{tot}}{\partial w} |}_{〈 w_{j} 〉} + {\frac{\partial^{2} Φ_{tot}}{\partial w^{2}} |}_{〈 w_{j} 〉}} .$ $Mathematical equation: ${\kappa _j} = \sqrt {{3 \over {\left\langle {{w_j}} \right\rangle }}{{\left. {{{\partial {\Phi _{{\rm{tot}}}}} \over {\partial w}}} \right|}_{\left\langle {{w_j}} \right\rangle }} + {{\left. {{{{\partial ^2}{\Phi _{{\rm{tot}}}}} \over {\partial {w^2}}}} \right|}_{\left\langle {{w_j}} \right\rangle }}} .$$ (10)

Note that when evaluating this property for galaxy disks, we assumed co-rotation of the gaseous and stellar disks, resulting in only a single value of these quantities per annuli.

2.3.2 Toomre stability parameter

The Toomre stability criterion (Safronov 1960; Toomre 1964) evaluates whether a disk is stable against collapse considering gravity, pressure, and shear. For the gaseous and stellar disks, we determined $Q_{g, j} = \frac{κ_{j} σ_{D, g, j}}{π G Σ_{g, j}},$ $Mathematical equation: ${Q_{g,j}} = {{{\kappa _j}{\sigma _{D,g,j}}} \over {\pi G{\Sigma _{g,j}}}},$$ (11a) $Q_{s, j} = \frac{κ_{j} σ_{D, s, j}}{3.36 G Σ_{s, j}},$ $Mathematical equation: ${Q_{{\rm{s,}}j}} = {{{\kappa _j}{\sigma _{{\rm{D, s,}}j}}} \over {3.36G{\Sigma _{{\rm{s,}}j}}}},$$ (11b)

where Q > 1 for a stable disk. Here we introduced for the gaseous and stellar disks, respectively, the surface densities Σ_g / Σ_s and the velocity dispersions σ_D,g / σ_D,s.

We followed Bottema (1993, but see also van der Kruit & Freeman 2011) to calculate the velocity dispersion of the stars³ as $σ_{D, s, j} = \frac{v_{c, j}}{2} \exp (- \frac{〈 w_{j} 〉}{2 w_{D}}),$ $Mathematical equation: ${\sigma _{{\rm{D, s,}}j}} = {{{v_{c,j}}} \over 2}\exp \left( { - {{\left\langle {{w_j}} \right\rangle } \over {2{w_{\rm{D}}}}}} \right),$$ (12)

with circular velocity $v_{c, j} = \sqrt{〈 w_{j} 〉 {\frac{\partial Φ_{tot}}{\partial w} |}_{〈 w_{j} 〉}} .$ $Mathematical equation: ${v_{c,j}} = \sqrt {\left\langle {{w_j}} \right\rangle {{\left. {{{\partial {\Phi _{{\rm{tot}}}}} \over {\partial w}}} \right|}_{\left\langle {{w_j}} \right\rangle }}} .$$ (13)

We assumed that the velocity dispersion of the cold gas equals the speed of sound of the interstellar medium, which correlates with the star formation surface density, i.e. $c_{s, cold, j} = α_{cold} + β_{cold} {(\frac{Σ_{SFR, j}}{M_{⊙} {kpc}^{- 2} {yr}^{- 1}})}^{γ_{cold}},$ $Mathematical equation: ${c_{{\rm{s, cold,}}j}} = {\alpha _{{\rm{cold}}}} + {\beta _{{\rm{cold}}}}{\left( {{{{\Sigma _{{\rm{SFR,}}j}}} \over {{{\rm{M}}_ \odot }{\rm{kp}}{{\rm{c}}^{ - 2}}{\rm{y}}{{\rm{r}}^{ - 1}}}}} \right)^{{\gamma _{{\rm{cold}}}}}},$$ (14)

with free parameters α_cold, β_cold, and γ_cold. This relationship shows significant scatter in observations across various redshifts (see e.g. Lehnert et al. 2009; Genzel et al. 2011; Green et al. 2014; Krumholz & Burkhart 2016; Zhou et al. 2017; Krumholz et al. 2018; Mai et al. 2024, and references therein). In our fiducial model we assume the parameter values α_cold = 5 km s⁻¹, β_cold = 20 km s⁻¹, and γ_cold = 1/3, i.e. a turbulence-dominated energy dissipation prescription for the cold gas (Zhou et al. 2017). We explore the impact of different values in Appendix A.1.

It is well known that the gaseous and stellar disks interact dynamically (e.g. Lin & Shu 1966; Bertin & Romeo 1988) and that, for example, the stability of a stellar disk may be impacted by even small amounts of gas. To retain the same guidelines for the Toomre stability parameter in the prescription provided by Kruijssen (2012), we followed the approach by Romeo & Wiegert (2011) and computed an ‘effective’ Toomre parameter as $Q_{eff, j}^{- 1} = {\begin{array}{l} ψ_{Q, j} Q_{s, j}^{- 1} + Q_{g, j}^{- 1} if Q_{s, j} \geq Q_{g, j}, \\ Q_{s, j}^{- 1} + ψ_{Q, j} Q_{g, j}^{- 1} otherwise, \end{array}$ $Mathematical equation: $Q_{{\rm{eff,}}j}^{ - 1} = \left\{ {\matrix{{{\psi _{Q,j}}Q_{s,j}^{ - 1} + Q_{g,j}^{ - 1}\quad {\rm{if}}\quad {Q_{s,j}} \ge {Q_{g,j}},} \hfill \cr {Q_{s,j}^{ - 1} + {\psi _{Q,j}}Q_{g,j}^{ - 1}\quad {\rm{otherwise,}}} \hfill \cr } } \right.$$ (15)

with the weighting factor $ψ_{Q, j} = 2 \frac{σ_{D, s, j} σ_{D, g, j}}{σ_{D, s, j}^{2} + σ_{D, g, j}^{2}} .$ $Mathematical equation: ${\psi _{Q,j}} = 2{{{\sigma _{{\rm{D, s,}}j}}{\sigma _{{\rm{D, g,}}j}}} \over {\sigma _{{\rm{D, s,}}j}^2 + \sigma _{{\rm{D, g,}}j}^2}}.$$ (16)

In Figure 2 we show an overview of the cold gas surface mass density, the epicyclic frequency, and the Toomre stability parameter from running our model on tree-files of the Millennium simulation. For comparison, we added the position of the solar neighbourhood.

The epicyclic frequency and cold gas surface density decrease at larger galactocentric distances. The Toomre parameter shows a more complex behaviour, typically ranging between one and ten, except for the centre of elliptical galaxies and the largest distances where it increases to higher values. The former effect is caused by the importance of the bulge component which is weaker in spirals. At large radii the disks of low-mass galaxies are barely populated with gas and stars, which drive Q_eff to large values. For the same reason the Toomre parameter drops for more massive galaxies, which is reflected by an increase in the disk’s scale length. Overall we found good agreement between our simulated galaxies and the solar neighbourhood.

Fig.2

Epicyclic frequency, cold gas surface mass density, and the Toomre stability parameter as a function of galactocentric distance for disk- (blue) and bulge-dominated (red) galaxies, defined as having a bulge-to-total stellar mass ratio of B/T < 0.2 and B/T ≥ 0.9, respectively. Vertical lines show the annuli’s outer radii, as defined in Eq. (1). We add for comparison the value of the solar neighbourhood: We calculated κ_D,⊙ ≈ 0.046 Myr⁻¹, as derived from the Oort constants A = 15.6 km s⁻¹ kpc⁻¹ and B = −15.8 km s⁻¹ kpc⁻¹ taken from Guo & Qi (2023); Σ_g,⊙ ≈ 13 M_⊙ pc⁻² from Flynn et al. (2006); and Q_eff,⊙ ≈ 1.7 (Binney & Tremaine 2008, with Q_s,⊙ ≈ 2.7 and Q_g,⊙ ≈ 1.5), a typical value for disks (e.g. Rafikov 2001; Leroy et al. 2008; Feng et al. 2014; Westfall et al. 2014). Note that we did not calculate the Toomre stability parameter for annuli, where its surface density drops below 1 M_⊙ pc⁻² as we do not expect the formation of any star clusters at such low gaseous densities.

2.3.3 Bound fraction of star formation

The bound fraction of newly formed stars in star clusters is closely related to the cluster formation efficiency (Bastian 2008; Goddard et al. 2010; Adamo et al. 2011; Silva-Villa & Larsen 2011), which takes into account its survival rate during the first few million years. Although important in many aspects, it is unclear how the bound fraction and cluster formation efficiency are related to the interstellar medium and star formation (see Andersson et al. 2024, for a recent discussion on how feedback influences the cluster formation efficiency).

As mentioned above, we follow the model outlined by Kruijssen (2012) to estimate the bound fraction based on its epicyclic frequency, cold gas surface density, and Toomre stability parameter, all equated at the log-mean galactocentric distance. We briefly highlight key aspects of the model and refer the interested reader to the original work for a more detailed description.

Note that we did not directly determine the bound fraction during the execution of the simulation. Instead, to reduce the computational cost, we created lookup tables for a set of {Σ_g, Q_eff, κ}. In total, we utilised 30 lookup tables with varying values of Q_eff (between 0.5 and 100) and 500 values of Σ_g and κ each, resulting in 7.5 × 10⁶ data points. The simulation then determined the closest match in all three parameters and extracted the bound fraction from the table.

Following an extensive literature (e.g. Padoan et al. 1997; Vázquez-Semadeni et al. 1998; Ostriker et al. 2001; Kritsuk et al. 2007; Padoan & Nordlund 2011; Kritsuk et al. 2017; Burkhart 2018) we assumed that the density contrast of the interstellar medium follows a log-normal distribution of the form $d p_{j} = \frac{1}{\sqrt{2 π} ς_{ρ, j}} \exp [- {(\frac{\ln δ_{j} - \bar{\ln δ_{j}}}{\sqrt{2} ς_{ρ, j}})}^{2}] d (\ln δ_{j})$ $Mathematical equation: ${\rm{d}}{p_j} = {1 \over {\sqrt {2\pi } {\varsigma _{\rho ,j}}}}\exp \left[ { - {{\left( {{{\ln {\delta _j} - \overline {\ln {\delta _j}} } \over {\sqrt 2 {\varsigma _{\rho ,j}}}}} \right)}^2}} \right]{\rm{d}}(\ln {\delta _j})$$ (17)

with relative density δ_j = ρ/ρ_{ISM, j}, ***(21)*** $\bar{\ln δ_{j}} = - 0.5 ς_{ρ, j}^{2}$ $Mathematical equation: $\overline {\ln {\rm{ }}{\delta _j}} = - 0.5\varsigma _{\rho ,j}^2$$ (e.g. Vazquez-Semadeni 1994) and standard deviation $ς_{ρ, j} = \sqrt{\ln (1 + 3 γ_{ρ}^{2} ℳ_{cold, j}^{2})},$ $Mathematical equation: ${\varsigma _{\rho ,j}} = \sqrt {\ln (1 + 3\gamma _\rho ^2{\cal M}_{{\rm{cold,}}j}^2)} $$ (18)

where γ_ρ ≈ 0.5 (Nordlund & Padoan 1999; Padoan & Nordlund 2002). The Mach number of the cold gas is related to the cold gas surface density, the epicyclic frequency, and the Toomre stability parameter via $Μ_{cold, j} = \sqrt{2} ϕ_{\bar{P}, j}^{1 / 8} \frac{Q_{eff,} j Σ g, j}{κ j},$ $Mathematical equation: ${{\rm M}_{{\rm{cold, }}j}} = \sqrt 2 \phi _{{\rm{\bar P,}}\,j}^{1/8}{{{Q_{{\rm{eff,}}}}\,j\Sigma {\rm{g,}}\,j} \over {\kappa j}},$$ (19)

where ***(24)*** $ϕ_{\bar{P}, j}$ $Mathematical equation: ${\phi _{{\rm{\bar P}},j}}$$ is the ratio of the mean pressure of a gaseous cloud related to the pressure at its surface (Krumholz & McKee 2005), with typical values close to two in dense regions (Heyer et al. 2004; Rosolowsky & Blitz 2005; Schuster et al. 2007; Colombo et al. 2014). Notice that the parameter ***(25)*** $ϕ_{\bar{P}, j}$ $Mathematical equation: ${\phi _{{\rm{\bar P}},j}}$$ is directly related to the fraction of cold gas contained in GMCs as ***(26)*** $ϕ_{\bar{P}, j} \approx 10 - 8 \times f_{GMC, j}$ $Mathematical equation: ${\phi _{{\rm{\bar P}},j}} \approx 10 - 8 \times {f_{{\rm{GMC}},j}}$$ . We assumed that this fraction only depends on the cold gas surface density (Krumholz & McKee 2005), i.e. $f_{GMC, j} = {[1 + 250 / {(Σ_{g, j} [M_{⊙} {pc}^{- 2}])}^{2}]}^{- 1} .$ $Mathematical equation: ${f_{{\rm{GMC,}}j}} = {\left[ {1 + 250/{{\left( {{\Sigma _{g,j}}[{{\rm{M}}_ \odot }{\rm{p}}{{\rm{c}}^{ - 2}}]} \right)}^2}} \right]^{ - 1}}.$$ (20)

Next, we needed to evaluate the minimum-value star formation efficiency. If star formation occurs on the free-fall timescale, this efficiency can be expressed as a combination of the specific SFR (sSFR_ff) and the ratio of the feedback timescale (t_fb) to the free-fall timescale (t_ff), i.e. $ε_{ff, j} = {sSFR}_{ff, j} t_{fb, j} / t_{ff, j} .$ $Mathematical equation: ${\varepsilon _{{\rm{ff,}}j}} = {\rm{sSF}}{{\rm{R}}_{{\rm{ff,}}j}}{t_{{\rm{fb,}}j}}/{t_{{\rm{ff,}}j}}.$$ (21)

For a spherical symmetric and homogeneous mass distribution, $t_{ff, j} = \sqrt{\frac{3 π}{32 G ρ_{ISM, j}}},$ $Mathematical equation: ${t_{{\rm{ff,}}j}} = \sqrt {{{3\pi } \over {32G{\rho _{{\rm{ISM,}}j}}}}} ,$$ (22)

where we determine the density of the interstellar medium as $ρ_{ISM, j} = \frac{ϕ_{P}}{2 π G} {(\frac{κ_{j}}{Q_{eff, j}})}^{2},$ $Mathematical equation: ${\rho _{{\rm{ISM,}}j}} = {{{\phi _{\rm{P}}}} \over {2\pi G}}{\left( {{{{\kappa _j}} \over {{Q_{{\rm{eff,}}j}}}}} \right)^2},$$ (23)

with ϕ_P ≈ 3 (see Appendix A of Krumholz & McKee 2005).

For the feedback timescale (Kruijssen 2012), $t_{fb, j} = \frac{t_{SN}}{2} [1 + \sqrt{1 + \frac{4 π^{2} G^{2} t_{ff, j}}{ϕ_{fb} {sSFR}_{ff, j} t_{SN}^{2}} {(\frac{Q_{eff, j} Σ_{g, j}}{κ_{j}})}^{2}}],$ $Mathematical equation: ${t_{{\rm{fb,}}j}} = {{{t_{{\rm{SN}}}}} \over 2}\left[ {1 + \sqrt {1 + {{4{\pi ^2}{G^2}{t_{{\rm{ff,}}j}}} \over {{\phi _{{\rm{fb}}}}{\rm{sSF}}{{\rm{R}}_{{\rm{ff,}}j}}t_{{\rm{SN}}}^2}}{{\left( {{{{Q_{{\rm{eff,}}j}}{\Sigma _{g,j}}} \over {{\kappa _j}}}} \right)}^2}} } \right],$$ (24)

where t_SN is the timescale for the first supernovae, which we assume to be 3 Myr, ϕ_fb = 5.28×10² pc²Myr⁻³ (Kruijssen 2012), and $\frac{{sSFR}_{ff, j}}{0.13} = 1 + erf [\frac{ς_{ρ, j}^{2} - \ln (0.68 α_{vir}^{2} ℳ_{cold, j}^{4})}{2^{3 / 2} ς_{ρ, j}}],$ $Mathematical equation: ${{{\rm{sSF}}{{\rm{R}}_{{\rm{ff,}}j}}} \over {0.13}} = 1 + {\rm{erf}}\left[ {{{\varsigma _{\rho ,j}^2 - \ln (0.68\alpha _{{\rm{vir}}}^2{\cal M}_{{\rm{cold,}}j}^4)} \over {{2^{3/2}}{\varsigma _{\rho ,j}}}}} \right],$$ (25)

with the virial parameter of GMCs α_vir (Larson 1981). It relates the mass, size, and velocity dispersion of a GMC and takes values between 10⁻¹ and 10¹ (e.g. Myers & Goodman 1988; Bertoldi & McKee 1992; Heyer et al. 2009; Dobbs et al. 2011; Hopkins et al. 2012). We set α_vir = 1.3 as proposed by McKee & Tan (2003).

In case star formation is less efficient than is assumed in Eq. (21) we take ϵ_{inc, j} = ϵ_{ff, j} × t_{ff, j}/10 Myr. If star formation is more efficient we set an upper bound of ϵ_{inc, j} = ϵ_max = 0.5 (Matzner & McKee 2000). The resulting effective star formation efficiency is the minimum of the above values, $ε_{j} = \min (ε_{ff, j}, ε_{inc, j}, ε_{max}) .$ $Mathematical equation: ${\varepsilon _j} = \min ({\varepsilon _{{\rm{ff,}}j}},{\varepsilon _{{\rm{inc,}}j}},{\varepsilon _{{\rm{max}}}}).$$ (26)

Finally, the bound fraction of star formation can be computed by the normalised integral of the probability density function of the contrast of the interstellar medium combined with the minimum star formation efficiency, $f_{bound, j} = ε_{max}^{- 1} \frac{\int_{- \infty}^{\infty} d δ_{j} ε_{j}^{2} (δ_{j}) δ_{j} (d p_{j} / d δ_{j})}{\int_{- \infty}^{\infty} d δ_{j} ε_{j (} δ_{j}) δ_{j} (d p_{j} / d δ_{j})} .$ $Mathematical equation: ${f_{{\rm{bound,}}j}} = \varepsilon _{{\rm{max}}}^{ - 1}{{\mathop \smallint \nolimits_{ - \infty }^\infty {\rm{d}}{\delta _j}\varepsilon _j^2({\delta _j}){\delta _j}({\rm{d}}{p_j}/{\rm{d}}{\delta _j})} \over {\mathop \smallint \nolimits_{ - \infty }^\infty {\rm{d}}{\delta _j}{\varepsilon _{j(}}{\delta _j}){\delta _j}({\rm{d}}{p_j}/{\rm{d}}{\delta _j})}}.$$ (27)

Figure 3 shows the bound fraction for different environments with Q_eff = 0.5 for the background. The z = 0 distribution of annuli from running L-Galaxies on halos identified in the Millennium and Millennium-II simulations is added on top and reveals a large range bound star formation that depends on the location within the galaxy: the innermost regions feature high surface densities and epicyclic frequencies (c.f. Figure 2) and have high bound fractions approaching one. In contrast, the outermost regions of galaxies have low surface densities and epicyclic frequencies, prohibiting the formation of bound structures. As a consequence, this result already predicts that massive star clusters at distances a few times the disk’s scale length likely originate from accreted satellite galaxies, assuming that heating processes within their host galaxy are insignificant.

While not shown in Figure 3, we find a decrease in bound fraction for an increasing Toomre stability parameter value when keeping the epicyclic frequency and cold gas surface density constant. This is related to the specific SFR, sSFR_ff, which decreases for an increasing Toomre parameter as per Eqs. (19) and (25), and because the gas disk becomes more stable with an increase in Q_eff. As a result the mid-plane density of the interstellar medium decreases, which, in turn, decreases the star formation efficiency and, thus, the bound fraction.

Fig. 3

Bound fraction of star formation, evaluated for Q_eff = 1.0, as a function of angular frequency and cold gas surface density. Solid blue and dashed grey contours give the smoothed distribution (with a standard deviation of 1 dex) of all annuli of all galaxies running L-Galaxies tree-files 0-9 and 40-79 on the Millennium and Millennium-II simulations, respectively. The location of the solar neighbourhood (see Figure 2 for details) is marked with a cross.

2.4 Initial star cluster properties

Here we detail the initial properties of star clusters that form in galaxy disks. Note that we deliberately exclude star clusters that form in the innermost ring with a radius of approximately 30 pc. These star clusters will either be disrupted quickly or merge to form a nuclear star cluster, which will be the subject of future work.

Another computational limit is the number of star clusters whose properties will be tracked over time. In this work, we focus on the most massive star clusters that survive until z = 0 and decide to completely ignore star clusters below 10³ M_⊙ (see below). In the simulation we track the properties of the 1000 most massive star clusters in the disk and halo, respectively, resulting in up to 2000 objects per galaxy. Although somewhat arbitrary, these numbers compare to the total number of star clusters of M 31 (potentially more than 1000; e.g. Barmby & Huchra 2001; Huxor et al. 2014) and is much larger than the number of known globular clusters for the Milky Way (about 200; see e.g. Minniti et al. 2017; Garro et al. 2024, and references therein).

2.4.1 Stellar mass

The mass of each star cluster is a random realisation of the underlying cluster initial mass function (CIMF), which we assumed follows a classical power-law distribution that is truncated at the upper mass end with an environmentally dependent mass-scale, i.e. $ξ_{j} (m_{c}) \propto m_{c}^{α_{CIMF}} \exp (- m_{c} / m_{cl,max, j}) .$ $Mathematical equation: ${\xi _j}({m_{\rm{c}}}) \propto m_{\rm{c}}^{{\alpha _{{\rm{CIMF}}}}}\exp ( - {m_{\rm{c}}}/{m_{{\rm{cl,max,}}j}}).$$ (28)

We assumed that α_CIMF = −2, motivated by observational studies of young star clusters in nearby galaxies (e.g. Zhang & Fall 1999; Bik et al. 2003; Hunter et al. 2003; McCrady & Graham 2007; Portegies Zwart et al. 2010; Emig et al. 2020; Levy et al. 2024). The truncation mass-scale is a product of the star formation efficiency ϵ_cloud, the bound fraction of star formation f_bound, the Toomre mass m_T, and the fraction of molecular gas that is critical to undergo gravitational collapse f_coll (Reina-Campos & Kruijssen 2017, but see also Kruijssen 2014 and Reina-Campos et al. 2022a), i.e. $m_{cl, max, j} = ε_{cloud} f_{bound, j} m_{T, j} f_{coll, j},$ $Mathematical equation: ${m_{{\rm{cl, max,}}j}} = {\varepsilon _{{\rm{cloud}}}}{f_{{\rm{bound,}}j}}{m_{{\rm{T,}}j}}{f_{{\rm{coll,}}j}},$$ (29)

and it results in a typical value of the order of 10⁵ M_⊙ at lower redshifts and up to 10⁹ M_⊙ in extreme cases. We assumed $m_{T, j} = 4 π^{5} G^{2} Σ_{g, j}^{3} Q_{eff, j}^{4} / κ_{j}^{4},$ $Mathematical equation: ${m_{{\rm{T,}}j}} = 4{\pi ^5}{G^2}\Sigma _{g,j}^3Q_{{\rm{eff,}}j}^4/\kappa _j^4,$$ (30)

and $f_{coll, j} = \min (^{1, t_{fb, j} / t_{ff, 2D, j}) 4},$ $Mathematical equation: ${f_{{\rm{coll,}}j}} = \min {(1,{t_{{\rm{fb,}}j}}/{t_{{\rm{ff, 2D,}}j}})^4},$$ (31)

for the Toomre mass and collapse fraction, respectively. The two-dimensional free-fall timescale is ***(39)*** $t_{ff,2D, j} = \sqrt{2 π} / κ_{j}$ $Mathematical equation: ${t_{{\rm{ff,2D, }}j}} = \sqrt {2\pi } /{\kappa _j}$$ .

In the above equations we assume a constant value of ϵ_cloud = 0.1 for star formation within a GMC, motivated by numerical results (e.g. Oklopčić et al. 2017; Chevance et al. 2020). Notice that this value is potentially smaller than the assumed star formation efficiency in the determination of the bound fraction in Eq. (26). A higher efficiency would result in a higher upper truncation mass-scale of the CIMF and could result in the formation of more massive star clusters. However, as we show in Appendix A.3 the main results of our work do not significantly change when assuming no upper truncation mass-scale but a pure power-law function instead.

For computational efficiency of the code, we consider only star clusters with initial masses in the range of 10⁴–10⁸ M_⊙. At the same time, we assume that the minimum star cluster mass that could, in principle, form is 10² M_⊙, meaning that we do not fully sample the CIMF. Therefore, we only sampled a total star cluster mass of f_sample f_bound SFR δt where $f_{sample} = \frac{\int_{10^{4} M_{⊙}}^{10^{8} M_{⊙}} m_{c} ξ (m_{c}) d m_{c}}{\int_{10^{2} M_{⊙}}^{10^{8} M_{⊙}} m_{c} ξ (m_{c}) d m_{c}} .$ $Mathematical equation: ${f_{{\rm{sample}}}} = {{\mathop \smallint \nolimits_{{{10}^4}\,{{\rm{M}}_ \odot }}^{{{10}^8}\,{{\rm{M}}_ \odot }} {m_{\rm{c}}}\xi ({m_{\rm{c}}}){\mkern 1mu} {\rm{d}}{m_{\rm{c}}}} \over {\mathop \smallint \nolimits_{{{10}^2}\,{{\rm{M}}_ \odot }}^{{{10}^8}\,{{\rm{M}}_ \odot }} {m_{\rm{c}}}\xi ({m_{\rm{c}}}){\mkern 1mu} {\rm{d}}{m_{\rm{c}}}}}.$$ (32)

In the case of a power-law CIMF with index −2, this ratio is two-thirds, i.e. one-third of the total mass in star clusters is contained in objects below 10⁴ M_⊙. For the more complex case of Eq. (28), we pre-compute the integral numerically for 21 equally spaced values between 10³ and 10⁹ M_⊙, resulting in, for example, f_sample ≈ 2.3 × 10⁻⁶ for m_cl,max = 10³ M_⊙. Afterwards, we randomly sample the CIMF with m_cl,max that agrees best with the computed value of Eq. (29).

Finally, to reduce computational cost we utilised lookup tables for initial star cluster masses. For simplicity, we generated 21 lookup tables for different values of m_cl,max equally spaced between 10³ and 10⁹ M_⊙ and with 10⁷ data points each.

2.4.2 Initial galactocentric distances

Within each annulus, we assume a uniform distribution of the initial galactocentric distances of newly formed star clusters. In our model, this assumption neglects both other parameters of the star cluster (like mass or radius), as well as their local environments, such as the density of cold gas in the disk.

2.4.3 Half-mass radius

The physical processes that govern the distribution of the initial half-mass radius of star clusters is still unknown, and many theoretically motivated and observational-based prescriptions seem to fail to reproduce the distribution at z = 0 in nearby galaxies (e.g. Reina-Campos et al. 2023a). For that reason we adopt a simplified prescription by using a constant initial value of r_c,h = 1.0 pc for all clusters, independent of mass and redshift of formation. We explore the impact of a different value and more complex prescriptions in Appendix A.2.

2.4.4 Tidal radius

The half-mass sizes of star clusters typically increase over time (c.f. Section 2.5.2) and are limited to the tidal radius where the gravitational acceleration of the cluster equals the tidal acceleration. The tidal field is directly related to the local gravitational potential and the tidal radius is related to the first eigenvector of the diagonalised tidal tensor. Assuming circular orbits and a mass concentration in the galaxy centre, we applied the definition of King (1962, but see also Renaud et al. 2011; Renaud 2018) and calculated $r_{c, t} = {[\frac{G m_{c}}{{(Ω^{2} - \partial^{2} Φ_{tot} / \partial w^{2})}_{w_{c}}}]}^{1 / 3},$ $Mathematical equation: ${r_{{\rm{c,}}t}} = {\left[ {{{G{m_{\rm{c}}}} \over {{{({\Omega ^2} - {\partial ^2}{\Phi _{{\rm{tot}}}}/\partial {w^2})}_{{w_{\rm{c}}}}}}}} \right]^{1/3}},$$ (33)

where Ω equals the angular frequency equated at the galactocentric distance of the star cluster. More generally, for each ring j, $Ω_{j} = \sqrt{\frac{1}{〈 w_{j} 〉} {\frac{\partial Φ_{tot}}{\partial w} |}_{〈 w_{j} 〉}} .$ $Mathematical equation: ${\Omega _j} = \sqrt {{1 \over {\left\langle {{w_j}} \right\rangle }}{{\left. {{{\partial {\Phi _{{\rm{tot}}}}} \over {\partial w}}} \right|}_{\left\langle {{w_j}} \right\rangle }}} .$$ (34)

Note that this prescription ignores the tidal effect of nearby baryonic over-densities, such as star forming regions or clouds, which may be dominant over the global galactic field. Other works that use hydrodynamical approaches (such as Reina-Campos et al. 2022a) determine a local tidal tensor from neighbouring cells, which is not possible in our model. We discuss this issue in Section 3.6.

2.4.5 Metallicity

For the metallicity (and individual elemental abundances of H, He, C, N, O, Ne, Mg, Si, S, Ca, Fe) of a star cluster we assume that it equals the metallicity of the cold gas in the ring it forms in, Z_c = Z_g, and that this value remains constant over the star cluster’s lifetime. The first assumption neglects any azimuthal variations of metallicity profiles which are known to exist in some galaxies due to asymmetric structures such as bars or spiral patters in the Milky Way (e.g. Poggio et al. 2022; Spina et al. 2022; Filion et al. 2023; Hawkins 2023; Hackshaw et al. 2024) and other spiral galaxies (e.g. Sánchez-Menguiano et al. 2016; Ho et al. 2017; Sánchez-Menguiano et al. 2018; Ho et al. 2019; Hwang et al. 2019; Metha et al. 2021, 2024) but see Kreckel et al. (2019) for counterexamples. Nevertheless, despite lacking such an implementation we argue in Section 3.5 that the metallicity distributions of our cluster populations are reasonable with respect to observational constraints.

2.5 Star cluster evolution

We detail in the next subsections the main processes that affect the evolution of star clusters. These include different mass loss rates, the expansion of the cluster’s half-mass radius, and prescription for a cluster’s relocation due to dynamical friction and galaxy mergers.

2.5.1 Mass loss

We considered three mechanisms for cluster mass loss: stellar evolution, tidal stripping due to an expanding cluster, and tidal shocks from interactions with GMCs, i.e. $\frac{d m_{c}}{d t} = {\frac{d m_{c}}{d t} |}_{ev} + {\frac{d m_{c}}{d t} |}_{rlx} + {\frac{d m_{c}}{d t} |}_{sh} .$ $Mathematical equation: ${{{\rm{d}}{m_{\rm{c}}}} \over {{\rm{d}}t}} = {\left. {{{{\rm{d}}{m_{\rm{c}}}} \over {{\rm{d}}t}}} \right|_{{\rm{ev}}}} + {\left. {{{{\rm{d}}{m_{\rm{c}}}} \over {{\rm{d}}t}}} \right|_{{\rm{rlx}}}} + {\left. {{{{\rm{d}}{m_{\rm{c}}}} \over {{\rm{d}}t}}} \right|_{{\rm{sh}}}}.$$ (35)

In the following paragraphs, we present how we evaluate each term.

Stellar evolution. To take into account cluster mass loss from stellar evolution we assume that cluster members represent a random realisation of a single stellar population assuming a Chabrier IMF (Chabrier 2003), resulting in varying expected lifetimes. Then, for an individual 10⁵ M_⊙ star cluster we utilise the “Stochastically Lighting Up Galaxies” library (da Silva et al. 2012, 2014; Krumholz et al. 2015) with non-rotating ‘Geneva’ 2013 stellar tracks (Schaller et al. 1992; Meynet et al. 1994; Ekström et al. 2012; Georgy et al. 2013) and fit a linear relationship to the retained mass as a function of time. The resulting mass-loss rate at time t in [yr] is ${\frac{d m_{c}}{d t} |}_{ev} = - \frac{0.13}{\ln (10)} \frac{m_{c}^{init}}{t}, for t \geq t_{SN} .$ $Mathematical equation: ${\left. {{{{\rm{d}}{m_{\rm{c}}}} \over {{\rm{d}}t}}} \right|_{{\rm{ev}}}} = - {{0.13} \over {\ln (10)}}{{m_{\rm{c}}^{{\rm{init}}}} \over t},\quad {\rm{for}}t \ge {t_{{\rm{SN}}}}.$$ (36)

We find that this relation holds, to the first order, irrespective of cluster metallicity. Note that the above determination does not take into account other mass-loss channels that we introduce below, which is why we updated the ‘initial’ star cluster mass, denoted here with ***(45)*** $m_{c}^{init}$ $Mathematical equation: $m_{\rm{c}}^{{\rm{init}}}$$ , at each time step.

Relaxation. Multi-body encounters between stars within a cluster result in energy transfer between the individual bodies and can cause stars to either orbit the star cluster’s centre at larger radii or leave the cluster completely in case its velocity exceeds the escape velocity. For bound stars, if the new orbit crosses the tidal radius, the star can be stripped from the cluster, resulting in an effective mass loss.

We considered this effect by following the literature (e.g. Spitzer 1940; Hénon 1961; Spitzer 1987; Lamers et al. 2005) and set ${\frac{d m_{c}}{d t} |}_{rlx} = - ξ_{rlx} \frac{m_{c}}{τ_{rlx}},$ $Mathematical equation: ${\left. {{{{\rm{d}}{m_{\rm{c}}}} \over {{\rm{d}}t}}} \right|_{{\rm{rlx}}}} = - {\xi _{{\rm{rlx}}}}{{{m_{\rm{c}}}} \over {{\tau _{{\rm{rlx}}}}}},$$ (37)

with a relaxation timescale $τ_{rlx} = 0.138 \frac{\sqrt{N}}{\ln (γ_{rlx} N)} \sqrt{\frac{r_{c, h}^{3}}{G 〈 m_{⋆} 〉}} .$ $Mathematical equation: ${\tau _{{\rm{rlx}}}} = 0.138{{\sqrt N } \over {\ln ({\gamma _{{\rm{rlx}}}}N)}}\sqrt {{{r_{{\rm{c,}}h}^3} \over {G\left\langle {{m_ \star }} \right\rangle }}} .$$ (38)

Here ⟨m_⋆⟩ is the average stellar mass of the star cluster, N = m_c / ⟨m_⋆⟩, and 0.07 ≲ γ_rlx ≲ 0.14 (Giersz & Heggie 1994). We assume a Chabrier (2003) initial stellar mass function, such that ⟨m_⋆⟩ = 0.42 M_⊙ and γ_rlx = 0.11. Finally, we choose ξ_rlx = 0.08 as suggested in the literature (e.g. Hénon 1961; Gieles et al. 2011; Gieles & Renaud 2016) for equal-mass cluster members, avoiding a proper treatment of the star cluster’s direct tidal environment (Alexander & Gieles 2012).

Tidal shocks. When a star cluster is located within the thin disk, i.e. has not been accreted during a galaxy merger event, it frequently interacts with GMCs if the fraction of cold gas bound within clouds is high. Depending on the impact parameter between the interaction the GMC can inject a significant amount of energy into the star cluster resulting in an increase in velocity dispersion and causing a fraction of the stars to escape the cluster as their velocity exceeds the cluster’s escape velocity.

To model this effect we approximate a cluster’s internal energy by assuming that it follows a Plummer-like density profile (Plummer 1911). Following the equations outlined in Kruijssen (2012), we set ${\frac{d m_{c}}{d t} |}_{sh} = - \frac{m_{c}}{τ_{sh}},$ $Mathematical equation: ${\left. {{{{\rm{d}}{m_{\rm{c}}}} \over {{\rm{d}}t}}} \right|_{{\rm{sh}}}} = - {{{m_{\rm{c}}}} \over {{\tau _{{\rm{sh}}}}}},$$ (39)

with $τ_{sh} = \frac{3 \sqrt{π}}{8 \sqrt{2^{2 / 3} - 1}} \frac{G}{g ϕ_{sh} ϕ_{ad} ϕ_{P} f_{sh} f_{GMC}} {[\frac{Q_{eff, j}}{κ_{j}}]}_{w_{c}}^{3} \frac{m_{c}}{r_{c, h}^{3}},$ $Mathematical equation: ${\tau _{{\rm{sh}}}} = {{3\sqrt \pi } \over {8\sqrt {{2^{2/3}} - 1} }}{G \over {g{\phi _{{\rm{sh}}}}{\phi _{{\rm{ad}}}}{\phi _{\rm{P}}}{f_{{\rm{sh}}}}{f_{{\rm{GMC}}}}}}\left[ {{{{Q_{{\rm{eff,}}j}}} \over {{\kappa _j}}}} \right]_{{w_{\rm{c}}}}^3{{{m_{\rm{c}}}} \over {r_{{\rm{c,}}h}^3}},$$ (40)

where g = 1.5, ϕ_sh = 2.8, f_sh = 3, and ϕ_ad = exp(−0.062). We used Eq. (20) to equate f_GMC for the annuli corresponding to the star cluster’s current galactocentric distance.

Note that only star clusters in galaxy disks are affected by tidal shocks due to encounters with GMCs, i.e. we set dm_c/dt|_sh = 0 for star clusters in a galaxy’s halo. Our implementation also does not couple the timescale of tidal shocks to the strength of the local tidal field tensor (see e.g. Alexander et al. 2014; Reina-Campos et al. 2022a, for details) because we do not attempt to model GMCs. However, because L-Galaxies already separates the atomic from the molecular gas in each annulus (see Henriques et al. 2020; Yates et al. 2021, but also Yates et al. 2024) future efforts may implement clouds and improve on the current prescription.

2.5.2 Radial expansion

Star clusters expand adiabatically due to contributions from mass loss from two-body interactions and tidal shocks (e.g. Reina-Campos et al. 2022a). This results in a radial expansion of the form $\frac{d r_{c, h}}{d t} = [(2 - f_{sh}^{- 1}) {\frac{d m_{c}}{d t} |}_{sh} + (2 - ζ ξ_{rlx}^{- 1}) {\frac{d m_{c}}{d t} |}_{rlx}] \frac{r_{c, h}}{m_{c}} .$ $Mathematical equation: ${{{\rm{d}}{r_{{\rm{c,}}h}}} \over {{\rm{d}}t}} = \left[ {(2 - f_{{\rm{sh}}}^{ - 1}){{\left. {{{{\rm{d}}{m_{\rm{c}}}} \over {{\rm{d}}t}}} \right|}_{{\rm{sh}}}} + (2 - \zeta \xi _{{\rm{rlx}}}^{ - 1}){{\left. {{{{\rm{d}}{m_{\rm{c}}}} \over {{\rm{d}}t}}} \right|}_{{\rm{rlx}}}}} \right]{{{r_{{\rm{c,}}h}}} \over {{m_{\rm{c}}}}}.$$ (41)

Following Gieles & Renaud (2016), we assumed that f_sh = 3, and following Gieles et al. (2011) $ζ ξ_{rlx}^{- 1} = \frac{5}{3} {[{(\frac{r_{c, h}}{r_{c, j}})}_{1} \frac{r_{c, j}}{r_{c, h}}]}^{3 / 2} \approx 0.092 {(\frac{r_{c, t}}{r_{c, h}})}^{3 / 2},$ $Mathematical equation: $\zeta \xi _{{\rm{rlx}}}^{ - 1} = {5 \over 3}{\left[ {{{\left( {{{{r_{{\rm{c,}}h}}} \over {{r_{{\rm{c,}}j}}}}} \right)}_1}{{{r_{{\rm{c,}}j}}} \over {{r_{{\rm{c,}}h}}}}} \right]^{3/2}} \approx 0.092{\left( {{{{r_{{\rm{c,}}t}}} \over {{r_{{\rm{c,}}h}}}}} \right)^{3/2}},$$ (42)

where we assumed in the approximation that the Jacobi radius equals the tidal radius and that (r_{c, h}/r_{c, j})₁ ≈ 0.145 (Hénon 1961). The particular choice of f_sh results in a contraction of a star cluster from encounters with GMCs as the mass loss from shocks is negative. A star cluster can only expand if the contribution of two-body relaxation dominates, which strictly requires ***(52)*** $ζ ξ_{rlx}^{- 1} > 2$ $Mathematical equation: $\zeta \xi _{{\rm{rlx}}}^{ - 1} > 2$$ or r_h ≲ 0.128 r_t. This condition is met in low-density environments, such as a galaxy’s halo, and we expect star clusters to contract in the inner parts of a galaxy’s disk.

2.5.3 Galactocentric migration

Without any external perturbations star clusters in a galaxy’s disk and halo migrate towards the inner regions through dynamical friction caused by interactions with constituents of three components: (1) field stars belonging to the stellar halo and bulge, (2) the gaseous halo, and (3) the dark matter halo. For the stellar and dark matter components, we follow the well-known Chandrasekhar (1943) prescription and assume an isotropic velocity distribution function of the components as well as circular orbits of the star clusters. As a consequence, a star cluster experiences a radial acceleration of $\frac{d v_{c}}{d t} = - 4 π G^{2} \frac{m_{c} ρ_{f}}{v_{c}^{2}} \ln Λ [erf (X) - \frac{2 X}{\sqrt{π}} \exp (- X^{2})],$ $Mathematical equation: ${{{\rm{d}}{v_{\rm{c}}}} \over {{\rm{d}}t}} = - 4\pi {{\rm{G}}^2}{{{m_{\rm{c}}}{\rho _{\rm{f}}}} \over {v_{\rm{c}}^2}}\ln \Lambda \left[ {{\rm{erf}}(X) - {{2X} \over {\sqrt \pi }}\exp ( - {X^2})} \right],$$ (43)

where X = v_c / (2σ_f) and ρ_f is the density of dark matter and stars in the halo evaluated at the position of the star cluster. The circular velocity of each object is computed via Eq. (13). For the Coulomb logarithm, we followed Binney & Tremaine (2008) and assumed $\ln Λ = \ln [\frac{w_{c}}{\max (r_{c, h}, G m_{c} / v_{typ}^{2})}],$ $Mathematical equation: $\ln \Lambda = \ln \left[ {{{{w_{\rm{c}}}} \over {\max ({r_{{\rm{c,}}h}},G{m_{\rm{c}}}/v_{{\rm{typ}}}^2)}}} \right],$$ (44)

where we replaced the maximum impact parameter with the galactocentric distance of the star clusters and assume for the typical velocity parameters of the dark matter halo, ***(55)*** $v_{ryp} = \sqrt{G M_{vir} / R_{vir}}$ $Mathematical equation: ${v_{{\rm{ryp}}}} = \sqrt {G{\rm{ }}{M_{{\rm{vir}}}}/{R_{{\rm{vir}}}}} $$ . For the velocity dispersion of halo stars we assume that their orbits are dominated by the underlying dark matter profile. Based on the Jeans equation for an isotropic system, Zentner & Bullock (2003) give the following approximation for the velocity dispersion, $σ_{H, s} = \frac{1.4393 x_{c}^{0.354}}{1 + 1.1756 x_{c}^{0.725}} v_{max},$ $Mathematical equation: ${\sigma _{{\rm{H, s}}}} = {{1.4393x_{\rm{c}}^{0.354}} \over {1 + 1.1756x_{\rm{c}}^{0.725}}}{v_{{\rm{max}}}},$$ (45)

with the maximum circular velocity within the dark matter halo v_max that we introduced in Eq. (3) and x_c = w_c/R_vir.

The radial acceleration caused by the hot gaseous component differs slightly from the one introduced in Eq. (43). Following Ostriker (1999) and Escala et al. (2004), we took $\frac{d v_{c}}{d t} = - \frac{4 π G^{2} m_{c} ρ_{H}}{v_{c}^{2}} g (ℳ_{hot}),$ $Mathematical equation: ${{{\rm{d}}{v_{\rm{c}}}} \over {{\rm{d}}t}} = - {{4\pi {{\rm{G}}^2}{m_{\rm{c}}}{\rho _{\rm{H}}}} \over {v_{\rm{c}}^2}}g({{\cal M}_{{\rm{hot}}}}),$$ (46)

where ρ_H is the density of the gaseous halo evaluated at the star clusters position and $g (ℳ_{hot}) = \frac{1}{2} \times {\begin{array}{l} \ln (\frac{1 + ℳ}{1 - ℳ}) - 2 ℳ if ℳ < 1.0, \\ \ln (1 - ℳ^{- 2}) + 2 \ln Λ if ℳ > 1.0. \end{array}$ $Mathematical equation: $g({{\cal M}_{{\rm{hot}}}}) = {1 \over 2} \times \left\{ {\matrix{{\ln \left( {{{1 + {\cal M}} \over {1 - {\cal M}}}} \right) - 2{\cal M}\quad {\rm{if}}{\cal M} < 1.0,} \hfill \cr {\ln (1 - {{\cal M}^{ - 2}}) + 2\ln \Lambda \quad {\rm{if}}{\cal M} > 1.0.} \hfill \cr } } \right.$$ (47)

The parameter ln Λ is the Coulomb-logarithm first introduced in Eq. (44). We calculated the Mach number of a cluster within the hot gas halo as ℳ_hot = v_c / c_s,hot, where (Tanaka & Haiman 2009; Choksi et al. 2017) $\frac{c_{s, hot}}{[{km s}^{- 1}]} = 1.8 \sqrt{1 + z} {(\frac{M_{vir}}{10^{7} M_{⊙}})}^{1 / 3} {(\frac{Ω_{M, 0} h^{2}}{0.14})}^{1 / 6} .$ $Mathematical equation: ${{{c_{{\rm{s, hot}}}}} \over {[{\rm{km}}{{\rm{s}}^{ - 1}}]}} = 1.8\sqrt {1 + z} {\left( {{{{M_{{\rm{vir}}}}} \over {{{10}^7}{\rm{}}{{\rm{M}}_ \odot }}}} \right)^{1/3}}{\left( {{{{\Omega _{{\rm{M, 0}}}}{h^2}} \over {0.14}}} \right)^{1/6}}.$$ (48)

2.5.4 Re-distribution during galaxy mergers

Explanations for the observed presence of star clusters in galaxy halos generally favour (a) interactions with massive perturbers (such as other star clusters or GMCs) that heat the object from its birth-place in a disk, (b) star cluster formation in galaxy outskirts during tidal interactions, or (c) a direct re-distribution during galaxy mergers. The lower tidal field and absence of tidal shocks may be an important ingredient in cluster survival as indicated by Baumgardt & Hilker (2018) who find that the Milky Way’s halo hosts all Galactic star clusters older than approximately five Gyr. Therefore, re-positioning of star clusters is an essential physical ingredient in simulating star cluster properties.

We consider two different scenarios based on the ratio of the baryonic masses of the two galaxies that participate in a merger. Here we follow the prescription by L-Galaxies and assume that a major merger occurs if the mass ratio exceeds q ≥ 0.1. In this case, we assume that the disks of both galaxies are destroyed and that all star clusters from both galaxies are contained within the halo of the successor galaxy. For minor galaxy mergers (q < 0.1) we assume that the star cluster population of the more massive galaxy remains unaffected and that all accreted star clusters migrate into the halo of their new host.

Generally speaking, the resulting baryonic distribution after galaxy mergers is sensitive to the initial conditions such as the galaxy mass ratio and their respective positions, orientations, and velocity vectors and magnitudes to each other. We therefore make simplifying assumption for the new galactocentric distances of accreted star clusters.

In L-Galaxies a galaxy starts to disrupt once the baryonic mass components equal or exceed the remaining dark matter mass. Once this condition is met, the model computes the following dynamical friction timescale (Binney & Tremaine 2008), $τ_{fric} = α_{fric} \frac{V_{200} d_{sat}^{2}}{G M_{sat} \ln Λ_{fric}},$ $Mathematical equation: ${\tau _{{\rm{fric}}}} = {\alpha _{{\rm{fric}}}}{{{V_{200}}d_{{\rm{sat}}}^2} \over {G{M_{{\rm{sat}}}}\ln {\Lambda _{{\rm{fric}}}}}},$$ (49)

where V₂₀₀ is the velocity at distance R₂₀₀ from the centre of the central halo, M_sat is the total (baryonic and dark matter) mass of the to-be accreted satellite, and d_sat is the distance from the satellite galaxy to the central’s centre. The parameter α_fric = 2.4 ensures that the bright-end of the z = 0 galaxy luminosity function is recovered (De Lucia & Blaizot 2007). Finally, ln Λ_fric is the Coloumb logarithm, set to ln Λ_fric = ln(1 + M₂₀₀/M_sat), where M₂₀₀ is the enclosed mass within R₂₀₀ of the central galaxy.

The model assumes that the satellite follows a circular trajectory such that its distance from the centre of the central halo is reduced by a factor of ∆t/τ_fric, where ∆t is the time since τ_fric was first computed, incremented according to the time step of the simulation. Once ∆t ≥ τ_fric, the satellite is assumed to be fully accreted.

In our approach we make use of these two timescales (τ_fric = 0 and τ_fric = ∆t) by tracking the distance of the satellite between both times⁴. At the same time, once the code computes τ_fric, we determine the radius of the galaxy from the condition that the baryonic mass of the satellite galaxy matches the dark matter mass, under the assumption that the density profile of the dark matter within this radius remains unchanged from the stripping of the outer parts. We found this radius (w_max) by solving numerically the following equation, $\frac{M_{Bary}}{M_{DM, 0}} = \frac{\ln (1 + c_{vir} w_{max} / R_{vir}) - c_{vir} / (R_{vir} / w_{max} + c_{vir})}{\ln (1 + c_{vir}) - c_{vir} / (1 + c_{vir})},$ $Mathematical equation: ${{{M_{{\rm{Bary}}}}} \over {{M_{{\rm{DM, 0}}}}}} = {{\ln (1 + {c_{{\rm{vir}}}}{w_{{\rm{max}}}}/{R_{{\rm{vir}}}}) - {c_{{\rm{vir}}}}/({R_{{\rm{vir}}}}/{w_{{\rm{max}}}} + {c_{{\rm{vir}}}})} \over {\ln (1 + {c_{{\rm{vir}}}}) - {c_{{\rm{vir}}}}/(1 + {c_{{\rm{vir}}}})}},$$ (50)

where M_Bary is the sum of the gaseous and stellar disks, the bulge, and central black hole, and M_DM,0 is the dark matter mass of the satellite prior to infall.

For the second distance, we assumed that the galaxy is completely disrupted once τ_fric = ∆t, where the central parts of the satellite get stripped, i.e. w_min = 0. We then related the internal distances w_min and w_max to the distances of the satellite during disruption d_min and d_max and the position of a star cluster in the satellite at position ***(62)*** $w_{c}^{old}$ $Mathematical equation: $w_{\rm{c}}^{{\rm{old}}}$$ to compute the new galactocentric distance ***(63)*** $w_{c}^{new}$ $Mathematical equation: $w_{\rm{c}}^{{\rm{new}}}$$ within the central halo via $w_{c}^{new} = {(\frac{w_{c}^{old}}{w_{max} - w_{min}})}^{α_{w}} (d_{max} - d_{min}) + d_{min},$ $Mathematical equation: $w_{\rm{c}}^{{\rm{new}}} = {\left( {{{w_{\rm{c}}^{{\rm{old}}}} \over {{w_{{\rm{max}}}} - {w_{{\rm{min}}}}}}} \right)^{{\alpha _w}}}({d_{{\rm{max}}}} - {d_{{\rm{min}}}}) + {d_{{\rm{min}}}},$$ (51)

where we assumed in our fiducial model that α_w = 4. Thus, this prescription assumes that the most distant but still bound star clusters relocate to a distance d_max and the most central ones to d_min. Star clusters with ***(65)*** $w_{c}^{old} > w_{\max}$ $Mathematical equation: $w_{\rm{c}}^{{\rm{old}}} > {w_{\max }}$$ are stripped from the satellite halo before the dynamical friction timescale is computed in the simulation. Due to the large distances involved (typically larger than 1 Mpc) we assume that these star clusters become unbound from the central galaxy. We do not track their subsequent evolution.

Note our assumption that all star clusters survive the tidal shock experienced during galaxy mergers, i.e. resulting in a “survival fraction” of unity (see Kruijssen & Cooper 2012; De Lucia et al. 2024, for a different approach). Furthermore, we do not add another mass-loss term for this scenario.

3 Results

In this work we focus on basic properties of the star cluster populations at z = 0. In the following sections, we refer to a star cluster as young (old) if its age is τ_c < 0.3 Gyr (τ_c ≥ 6 Gyr), unless specified otherwise. These age cuts are used throughout the analysis, and while they are somewhat arbitrary, we chose them to be able to better compare to both observational results and predictions by other simulations. Disk-dominated (‘spiral’) galaxies were assumed to be the ones with a bulge-to-total stellar mass ratio of B/T < 0.2, whereas bulge-dominated (‘elliptical’) systems have B/T ≥ 0.7⁵.

We focus our analysis on the output of running the model on Millennium (Springel et al. 2005) tree-files 0-9 (out of 512 total) that contain 118 558 galaxies at z = 0 and provide us with a representative sub-sample. The model performed similarly when running on Millennium-II (Boylan-Kolchin et al. 2009) tree-files 40-79 (a representative sub-sample; out of 512), going down to lower galaxy stellar masses. We provide a brief comparison in each subsection in case of differences.

3.1 M_V-SFR relationship

A first test for our model is to reproduce the empirical relationship between the absolute V-band magnitude of the brightest young star cluster (local quantity within a galaxy) versus the host galaxy’s SFR (global quantity) within nearby disk-dominated galaxies (e.g. Larsen 2002; Bastian 2008; Larsen 2010). Since this relation is not used as input for the simulation it serves as both a check of the models capabilities and a test to explore secondary correlations with other third quantities.

For the simulated data, as we store star cluster masses, we need to convert to the Johnson-Cousins V-band. We perform the conversion by using the Python version (Johnson et al. 2023) of the “Flexible Stellar Population Synthesis” code (Conroy et al. 2009, 2010; Conroy & Gunn 2010b,a) that uses as input the metallicity and age of a stellar population and yields an absolute magnitude in selected filterbands. Furthermore, for the computation, we assume that all star clusters are composed of a single stellar population that follows a Chabrier (2003) stellar initial mass function. To compare to observations we compile the data from Johnson et al. (2000); Larsen (2002); Rafelski & Zaritsky (2005); Bastian (2008); Annibali et al. (2009); Goddard et al. (2010); Adamo et al. (2011); Annibali et al. (2011); Pasquali et al. (2011); Silva-Villa & Larsen (2011); Cook et al. (2012); Ryon et al. (2014); Whitmore et al. (2014); Adamo et al. (2015); Lim & Lee (2015); Cook et al. (2023) for a local sample of galaxies. Additionally, we use data from Maschmann et al. (2024); Thilker et al. (2025), representing the results for PHANGS galaxies (see Lee et al. 2022, 2023, for details). For that dataset specifically, we use stellar mass estimates from Emsellem et al. (2022) and SFR estimates from Sun et al. (2022). Otherwise, galaxy stellar masses are taken from the 50 Mpc catalogue provided by Ohlson et al. (2023).

We show the resulting parameter space for young massive star clusters in disk-dominated galaxies evaluated at z = 0 in Figure 4. Our results show an increasing star cluster mass with increasing SFR, in qualitative agreement with the observations, although with a slightly steeper slope. To quantify the level of (dis-)agreement we perform a linear fit to the data sample and obtain uncertainties through 10 000 Monte Carlo iterations. The resulting slope values are ***(66)*** $α_{sim} = - {2.466}_{- 0.006}^{+ 0.006}$ $Mathematical equation: ${\alpha _{{\rm{sim}}}} = - 2.466_{ - 0.006}^{ + 0.006}$$ for the simulated data (N_gal = 66292) and ***(67)*** $α_{obs} = - {2.02}_{- 0.09}^{+ 0.09}$ $Mathematical equation: ${\alpha _{{\rm{obs}}}} = - 2.02_{ - 0.09}^{ + 0.09}$$ for the observations (N_gal = 103). Despite small differences in the slope, most of our star clusters are located in galaxies with SFRs of the order of 10⁻¹ M_⊙ yr⁻¹ and have absolute V-band magnitudes of M_V ≈ −11, which is in excellent agreement with observations.

When running the model on Millennium-II data we find a better agreement for the lowest-mass (but still the most massive and young) star clusters towards lower galaxy SFRs, which is because of the higher mass resolution offered by the simulation. We observe that the slope values of the relationship changes: applying a galaxy mass cut of 10^9.5 M_⊙ results in a slope of ***(68)*** $α_{sim,mrii} = - {2.02}_{- 0.03}^{+ 0.03}$ $Mathematical equation: ${\alpha _{{\rm{sim,mrii}}}} = - 2.02_{ - 0.03}^{ + 0.03}$$ . This value matches the results for the observational dataset.

We identify a secondary dependence on galaxy stellar mass in panel (B) where the slope of the relationship is steeper for massive galaxies: when constraining the galaxy sample to stellar masses above 10¹⁰ M_⊙ the slope steepens to α_sim ≈ −3.16 (N_gal = 14091)⁶. This trend is mainly related to the upper truncation mass-scale for the CIMF, which decreases due to a decreasing Toomre mass and bound fraction. Galaxies with stellar masses below 10^9.5 M_⊙ have comparable Q_D,g values to more massive galaxies but their Q_D,s values are larger because of their disk mass is dominated by gas and not stars and ***(69)*** $Q_{D, i} \propto Σ_{D, i}^{- 1}$ $Mathematical equation: ${Q_{{\rm{D,}}i}} \propto \Sigma _{{\rm{D}},i}^{ - 1}$$ . As a consequence, the slope value of these low-mass galaxies is α_sim ≈ −2.36 (N_gal = 43211).

For the cluster formation efficiency, which equals the bound fraction introduced in Eq. (27) and the survival rate of star clusters during the initial few megayears (see Kruijssen & Cooper 2012; Kruijssen 2012, for details on this “cruel cradle effect”), we find an increase with SFR (panel C). This is expected given the direct relationship between the SFR surface density and cold gas surface density (c.f. Eq. (2)) and because f_bound positively correlates with Σ_g (c.f. Figure 3). Our results are in excellent agreement with literature data.

In the observational studies that we mention in Section 3 the cluster formation efficiency is computed by comparing the total mass in young (τ_c ≤ 10 Myr) star clusters with the current SFR (e.g. Goddard et al. 2010; Ryon et al. 2014; Hollyhead et al. 2016). Broad-band photometry and fits to the inferred spectral energy distribution yield estimates of the star cluster parameters whereas Hα or infrared fluxes estimate the SFR.

We conclude that this relationship is sensitive to the resolution of the simulation. Taking into account biases in the mass range of the selected galaxy samples is important.

Fig. 4

Absolute V-band magnitude of the youngest and most massive star cluster versus the galaxy-averaged SFR. The galaxy sample is limited to disk-dominated galaxies that have a bulge-to-total stellar mass ratio of B/T < 0.2. We compare our results to various observations of nearby disk-dominated galaxies (see main text for details). For both the simulated data and the observations, we set an age cut of τ_c ≤ 0.3 Gyr on the star clusters. Panel A: full observational and simulated data samples. For the simulated data, we show the 1σ, 2σ, and 3σ intervals using blue solid lines. Panel B: same as in the first panel but with all data points colour-coded according to the host galaxy’s stellar mass. If no stellar mass estimate is available for observational data points, we show them with grey symbols. Panel C: same as the central panel but the data points are colour-coded by the cluster formation efficiency, which is a combination of the bound fraction of star formation and the ‘cruel cradle effect’ (Kruijssen & Cooper 2012; Kruijssen 2012) that takes the interaction of a proto-star cluster with its natal environment and nearby GMCs into account. Note that the two outliers, NGC 1705 and NGC 5238, are starburst galaxies and that their massive star clusters were previously classified as nuclear star clusters (Pechetti et al. 2020; Hoyer et al. 2021). Nuclear star clusters often exhibit complex formation histories (e.g. Spengler et al. 2017; Kacharov et al. 2018; Fahrion et al. 2021) and cannot easily be compared to our simulated star clusters.

3.2 Stellar masses

We show in Figure 5 the one-dimensional histogram of star cluster masses at z = 0. The galaxies are split into different morphologies and mass bins, and star clusters are separated by their age and location (disk versus halo). There exist two sharp truncations at mass values of 10³ and 10⁴ M_⊙, which correspond to the cut-off and lowest initial star cluster masses of the simulation. Both features disappear towards more massive galaxies as they contain more than 1000 star clusters in disks, thus increasing the lowest star cluster mass that we evolve⁷.

Young star clusters (τ_c ≤ 300 Myr) follow approximately a power-law distribution that is truncated towards the high-mass end, according to the initial star cluster mass function. Higher-mass galaxies can feature more massive star clusters, which is both a statistical effect from sampling the CIMF as well as higher upper truncation masses from a larger Toomre mass and bound fraction (Eq. (29)). Apparent differences in the histograms for young star clusters between galaxy morphologies exist mostly because of the technical limitation of only tracking a certain number of star cluster per galaxy. L-Galaxies assumes that star formation only occurs in galaxy disks, and both morphological types are expected to have similar disk properties, leading to the expectation of similar power-law distributions of star cluster masses (see Henriques et al. 2015, 2020, for details)⁸.

Looking at old (τ_c > 6 Gyr) star clusters we find broad Gaussian-like distributions. They broaden and shift towards higher star cluster masses due to a typically higher upper truncation mass of the initial star cluster mass function and the positive correlation between the number of star clusters and a galaxy’s mass. There exist differences between galaxy morphologies, which become most striking in the most massive galaxies. The masses of gaseous disks are elevated in disk-dominated systems compared to their bulge-dominated counterparts resulting in a higher number of massive star clusters. This observation is reversed in galaxy halos, which is related to the merger history of the two morphological types: galaxy mergers, both minor and major, increase the population of star clusters in a galaxy’s halo and this effect occurred more often in z = 0 bulge- than in disk-dominated galaxies.

Finally, we compare our simulated star clusters to observational constraints from Brown & Gnedin (2021) who analysed the structural properties of young star clusters in a set of nearby galaxies from the LEGUS programme (Calzetti et al. 2015). For the comparison, we apply the same galaxy mass and star cluster age constraints as for the simulated data. For the comparison with our simulated data, we split the literature values into the same galaxy mass bins (M_⋆/M_⊙ < 5 × 10⁹, 5 × 10⁹ ≤ M_⋆/M_⊙ < 10¹⁰, and 10¹⁰ ≤ M_⋆/M_⊙) and applied the same constraints on star cluster ages (τ_c ≤ 300 Myr for young and τ_c > 6 Gyr for old objects).

Star cluster ages in the analysis from Brown & Gnedin (2021) stem from fitting the spectral energy distribution; see details in Adamo et al. (2017) and Ryon et al. (2017)⁹. For dwarf galaxies, we find excellent agreement between the stellar mass distributions of the populations of simulated and observed star clusters down to the sampling limit of 10⁴ M_⊙. In more massive galaxies the comparison becomes more challenging due to a limited number of star cluster that we consider per galaxy. Nevertheless, we find that our simulated galaxies form more massive star clusters than observed in nearby systems, which could be related to multiple effects. First, to some degree, the comparison could suffer from poor number statistics from the observations, which does not populate well the high-mass end of star clusters. Second, our data contain highly star forming galaxies with SFRs above 10 M_⊙ yr⁻¹ at the high galaxy-mass end, which results in a larger probability to form massive star clusters. While our assumption for the formation of star clusters, such as the upper truncation mass of the CIMF, could be inaccurate and lead to elevated star cluster masses, the agreement in the M_V-SFR relationship that we presented in the previous section would indicate that we sample well the initial star cluster masses. Finally, the evolution of the young star clusters is sensitive to their initial properties and affected by their direct surroundings. As we present in the next section, there is a disagreement between the half-mass radii of the simulated and observed star clusters, which could drive some discrepancy between the observed mass functions, given the coupled differential equations of the evolution of a star cluster’s mass and radius.

Fig. 5

Estimated probability density function values of binned star cluster masses in disk- (B/T < 0.2; blue colour) and bulge-dominated (B/T > 0.7; red colour) galaxies at z = 0. We track up to 2000 star clusters per galaxy, split equally between star clusters located in the disk (top row) and halo (bottom row). Host galaxies are separated into three mass bins. Additionally, we separate between young (τ_c ≤ 300 Myr; solid lines) and old (τ_c > 6 Gyr; shaded regions) star clusters. We sample the CIMF with an initial mass of 10⁴ M_⊙ and discard star clusters with masses below 10³ M_⊙. The black solid line in the top left panel presents the asymptotic power-law behaviour towards m_c = 0. We compare the population of simulated star clusters to the observational resulted of Brown & Gnedin (2021, orange colour). Each column shows the same dataset as the location of the star clusters (disk versus halo) is not specified by the authors.

Fig. 6

Same as Figure 5 but for star cluster half-mass radii. The double-peaked structure of the half-mass radii in a galaxy’s halo is a result of the prescription for star cluster re-distribution during galaxy mergers (see Sections 2.5.4 and 3.4 below).

3.3 Half-mass radii

We present in Figure 6 the distribution of half-mass radii of star clusters at z = 0. The separation into different star cluster ages and galaxy masses follows the previous section.

We find that the half-mass distribution of young star clusters peak around the initial value of 1 pc and extends in both directions. Going from low- to high-mass galaxies we find that young star clusters become larger, which is a consequence of the increased average star cluster mass (see Eq. (41) and the related individual mass-loss terms). For the same reason, we observe the opposite trend with galaxy morphology between young star clusters in the disk and in the halo (compare to Figure 5)¹⁰.

Older star clusters in galaxy disks show a broader distribution of half-mass radii, which is related to both a broad distribution of their masses and a great variety of their environments, which affect the contributions of tidal shocks to the radial expansion or contraction (Eq. (40)). Additionally, the half-mass radii increase, ranging between a few parsecs in dwarf galaxies to ten parsecs in the most massive galaxies. The half-mass radii of star clusters in disks is limited due to the influence of tidal shocks on the mass: if a star cluster expands significantly, the mass-loss term for tidal shocks increases, thus reducing both the stellar mass and the leading to a contraction of the star cluster.

This coupling then leads to either the destruction of star clusters that experience repeated episodes of radial expansion and contraction or restrains the half-mass radii to small values.

Older star clusters in a galaxy halo, especially in bulge-dominated galaxies, show a peculiar distribution with a peak around a few parsecs and a separate extended tail towards large half-mass radii. Both trends are a direct consequence of our assumed prescription for the re-distribution of star clusters during galaxy mergers, which we detail in the next section. The peak around a few parsecs comes from major merger events and matches the star cluster distribution in disks, i.e. the population of star clusters, which were heated during the galaxy-galaxy interaction. In contrast, star clusters with large half-mass radii got accreted during minor merger events and are located at great distances from their new host galaxy. These star clusters do not suffer from tidal shocks. Instead, the large tidal radius of the star clusters enables two-body relaxation to increase their half-mass size to large values, often exceeding 100 pc (see Eqs. (41) and (42)).

We compare our simulated young star clusters to the results of Brown & Gnedin (2021) after applying the same constraints on the star cluster age and galaxy mass. The simulated star clusters show a much narrower distribution that does not extend to large radii. Instead, the half-mass radii of observed young star cluster fit well to the distribution of old star clusters in our simulation. If the data by Brown & Gnedin (2021) present a representative distribution of half-mass radii shortly after cluster birth then they indicate that the radii are potentially constrained by the clusters environment and may evolve quickly (see e.g. Banerjee & Kroupa 2017, for N-body simulations). Evidence for the latter comes from hydrodynamical simulations of Lahén et al. (2025) that show a rapid half-mass radius evolution from ≈0.1 to ≈1 pc over 100 Myr. We explore different prescriptions for initial half-mass radii in Appendix A.2 but note that Reina-Campos et al. (2023a) found that none of their prescriptions, including one that depends on the local environment, can recover the z = 0 distribution of the Milky Way.

Fig. 7

3.4 Galactocentric distances

We show the distribution of galactocentric distances of star clusters in Figure 7. The separation into different star cluster ages and galaxy masses follows the previous two sections.

Looking first at star clusters within a galaxy disk, we find the ‘step-like’ behaviour of young star clusters, which shows the assumed logarithmically flat distribution of star cluster in galactocentric distances. Deviations from this distribution only occur for young star clusters in the central few tens of parsecs, where the local density is highest, the most massive star clusters out to a few hundred parsecs, or old star clusters. The latter effect is clearly visible by considering only star clusters older than τ_c > 6 Gyr. Within ≈1000 pc star clusters efficiently sink towards the centre and either merge within the central few parsecs to form a nuclear star cluster, whose formation we do not consider here, or get tidally destroyed due to the increasing gas density and fraction of GMCs.

For star clusters in the halo we find two separate distributions. First, we find the same ‘step-like’ behaviour as for star clusters in the disk, which come from major mergers: as mentioned in Section 2.5.4 we assume that star clusters located in the disk of the more massive galaxy get heated to the halo but retain their galactocentric distance during major galaxy mergers. This effect, therefore, is mostly seen in bulge-dominated galaxies given that most z = 0 disk-dominated galaxies never experienced a major merger. In contrast, we assume that the less massive galaxy gets tidally disrupted in (minor and major) galaxy mergers, resulting in a large range of galactocentric distances.

Given our assumptions on how star clusters accrete, we find that the in situ star cluster population, that is present in the disk, dominates in number over the ex-situ star cluster population within a few kiloparsecs to ≈10 kpc, depending on the host galaxy mass and type. These values are in rough agreement with the values presented in Keller et al. (2020) for the E-MOSAICS simulation; however, the shape of the distribution of star cluster galactocentric distances is clearly different. A more realistic ‘smooth’ star cluster distribution as a function of distance that originates from major mergers requires a more accurate treatment for the more massive galaxy during galaxy mergers within L-Galaxies, which is beyond the scope of this work.

Fig. 8

Mean star cluster metallicity per galaxy versus host galaxy stellar mass, separated into disk- (left panel) and bulge-dominated (right panel) galaxies. Strong (faint) contours give the 1, 2, and 3σ distribution for old (young) star clusters. The dashed grey line gives the fiducial model of Pfeffer et al. (2023), which assumes, similar to our model, an environmentally dependent prescription for the upper truncation mass of the CIMF and the cluster formation efficiency. The solid black line gives the empirical relationship for ellipticals in the nearby Virgo galaxy cluster (Peng et al. 2006). Data for the Milky Way and M 31 stem from a self-compiled data table that will be presented in future work. Other data points come from Usher et al. (2012), Sesto et al. (2018), and Fahrion et al. (2020). The grey shaded area marks the ‘lower-limit floor’ at ⟨[Fe/H]⟩ = −2.5 of observed star cluster metallicities in other galaxies (e.g. Beasley et al. 2019).

3.5 Metallicity distributions

3.5.1 Mean cluster metallicity

We show in Figure 8 the mean metallicity of a galaxy’s star cluster population, as traced by their iron-abundances, versus galaxy stellar mass. The iron abundance values are calculated by the galactic chemical evolution model introduced into L-Galaxies by Yates et al. (2013) that takes into account contributions from stellar winds and (type-Ia and -II) supernovae with different delay-time-distributions. We split the galaxies by their morphological type (disk-dominated versus elliptical; see above) and separate the star cluster population by their age (young versus old; see above).

Irrespective of the star clusters age and galaxy morphology we find that the mean metallicity increases as a function of galaxy mass, similar to the mass-metallicity relationship for galaxies (as traced via oxygen-abundances; e.g. Tremonti et al. 2004; Kewley & Ellison 2008; Torrey et al. 2019; Sanders et al. 2021). Despite some overlap, younger star clusters have higher metallicity than their older counterparts at fixed stellar mass for both galaxy types. This difference appears to be more significant for disk- than for bulge-dominated galaxies, which is related to the origin of the cold gas that forms stars: in disk-dominated systems stars are predominantly formed from gas that has been continuously enriched by the above-mentioned stellar winds and supernovae channels resulting in metal-rich young star clusters. In bulge-dominated galaxies that form stars (and star clusters) a significant fraction of the cold gas comes from accreted lower-mass systems that contain relatively metal-poor gas.

Towards the low-mass end, we find that a significant fraction of galaxies (32.1 and 53.0% for disk- and bulge-dominated galaxies below 10¹⁰ M_⊙, respectively) have at least one star cluster with a metallicity value ⟨[Fe/H]⟩ ≤ −2.5. In contrast, no disk- and only 0.003% of bulge-dominated galaxies have a mean star cluster metallicity below −2.5. This threshold value of −2.5 is often used to indicate a “metallicity floor” due to an apparent lack of globular clusters in nearby galaxies below this value (Beasley et al. 2019)¹¹. The contours of our simulated data drop down to extremely low metallicities at the mass-resolution limit. When running the model on the Millennium-II dataset, the star cluster distributions follow more closely the scaling relation by Peng et al. (2006), indicating that resolution effects play a role at the low-mass end.

We compare our results to observational data. For disk-dominated galaxies, we consider the Milky Way and M 31 as they have the most robust and quantitative measurements of globular cluster metallicities. To obtain the mean metallicity values, we collected the cluster information from a diverse set of literature that we present in a future paper (another table is presented in Pace 2025). Here we find that our simulated disk-dominated galaxies have, on average, about 0.3 dex higher median metallicity than the observed values in the Milky Way, which might be related to a (1) a lack of accreted dwarf galaxies that would contribute mainly low-metallicity star clusters, (2) a too aggressive prescription for the destruction of star clusters with low metallicity or old age, or (3) a computational limitation in that we only store the properties of the 1000 most massive star clusters in a disk. Pfeffer et al. (2023) find a better agreement with observational constraints when limiting their sample to a galactocentric distance of 20 kpc, thus excluding preferentially low-metallicity star clusters. We explore this issue again in a future paper.

For bulge-dominated galaxies we collect data from Usher et al. (2012); Sesto et al. (2018); Fahrion et al. (2020) and take the empirical scaling relation from Peng et al. (2006). Our simulated star cluster populations show excellent agreement with the observations across the whole mass scale. A couple of literature data points lie outside the 3σ contours, which could be related to either poor number statistics of the or a bias in the ages of the star clusters in the observations with respect to our simulated galaxies.

3.5.2 Bimodality

An extensive set of literature work argues that the star cluster population of many, perhaps all, galaxies shows a bi- or multi-modal colour or metallicity distribution (e.g. Cohen et al. 1998; Kundu & Whitmore 2001; Beasley et al. 2008; Brodie et al. 2012; Blom et al. 2012; Usher et al. 2012, 2013; Escudero et al. 2015; Caldwell & Romanowsky 2016; Bassino & Caso 2017; Villaume et al. 2019; Fahrion et al. 2020; Hixenbaugh et al. 2022; Lomelí-Núñez et al. 2024). An often referred to explanation is that the bi-modality is a result of different origins of clusters with the bluer (more metal-poor) population having an ex-situ origin (e.g. Strader et al. 2005; Brodie & Strader 2006; Katz & Ricotti 2013; Tonini 2013). Numerical work has argued that the bi- or multi-modality could be a result of different epochs of star cluster formation; different cluster origins, as suggested by the above literature; or details related to cluster formation and destruction, or some combination thereof (see e.g. Kruijssen 2015; Choksi & Gnedin 2019; Chen et al. 2025). However, as pointed out by the observational work of Pastorello et al. (2015) and the numerical findings by Pfeffer et al. (2023) the bimodality may only be present in a minority (≲50%) of systems.

To determine the bimodality in our simulation we use an Bayesian Gaussian Mixture Model approach with Dirichlet initial conditions for all galaxies than contain more than 30 clusters.

Following Muratov & Gnedin (2010) and Pfeffer et al. (2023), we first determined $- 2 \ln λ = - 2 \ln [\frac{\max (ℒ_{1})}{\max (ℒ_{2})}],$ $Mathematical equation: $ - 2\ln \lambda = - 2\ln \left[ {{{\max ({{\cal L}_1})} \over {\max ({{\cal L}_2})}}} \right],$$ (52)

where max(ℒ_j) is the maximum value of the likelihood function evaluated over the metallicity distribution when considering j number of Gaussians (either one or two). Afterwards, we perform 100 bootstrap iterations to evaluate the probability that the solution is bi-modal with a probability threshold of 90%, as chosen in Pfeffer et al. (2023)¹². Finally, we calculated the weighted distance between the two Gaussians, $D_{G} = \frac{| μ_{1} - μ_{2} |}{\sqrt{σ_{1}^{2} + σ_{2}^{2}}} \overset{!}{\geq} \sqrt{2},$ $Mathematical equation: ${D_{\rm{G}}} = {{|{\mu _1} - {\mu _2}|} \over {\sqrt {\sigma _1^2 + \sigma _2^2} }}\mathop \ge \limits^! \sqrt 2 ,$$ (53)

where µ_j and σ_j are the mean and standard deviations of Gaussian j, respectively. If ***(72)*** $D_{G} < \sqrt{2}$ $Mathematical equation: ${D_{\rm{G}}} < \sqrt 2 $$ , we classify the distribution as unimodal, as adopted in Muratov & Gnedin (2010) and Pfeffer et al. (2023)¹³.

We show the bimodality of the star cluster distribution as a function of galaxy stellar mass in Figure 9 where we split all galaxies into disk- and bulge-dominated. Irrespective of galaxy morphology we find that the bimodal fraction (i.e. the fraction of galaxies that exhibits a bimodal metallicity distribution) is constrained to values ranging between ≈20 and ≈50%.

In disk-dominated galaxies, the bimodality of the star cluster population decreases with increasing galaxy stellar mass, which is attributed to both a continuous formation of star clusters of increasing metallicity, which does not show an intrinsic bimodal distribution at z = 0, and an increasing number of minor merger events, which bring in additional star cluster populations and diffuse any potential bimodality that would exist in the data. The latter argument can be used to explain why the bimodality decreases for the most massive bulge-dominated galaxies, where these systems experienced, on average, more minor and major mergers than their disk-dominated counterparts. As a result the overall metallicity distribution at z = 0 extends over a large range and individual sub-structures or features become diffused. This effect can be seen in observations as well, such as massive galaxies such as M 87 (see Fig. 20 in Cohen et al. 1998) or massive galaxy clusters (e.g. Harris et al. 2017).

For bulge-dominated systems specifically, we find a dip in the bimodal fraction around 8 × 10⁹ M_⊙ by up to ten percent. This decrease is related to an increased number of accreted star clusters in comparison to in situ star clusters, and is not a consequence of a sudden change in the number or type of galaxy mergers. Towards higher galaxy masses, the in situ fraction of star clusters increase again and result in an increase of the bimodality, continuing to a few times 10¹⁰ M_⊙ after which the number of accreted star clusters dominate again.

We compare our results to the E-MOSAICS simulation where Pfeffer et al. (2023) performed part of the metallicity analysis. Note that, for the comparison, we split the star clusters into two age bins, τ_c ≥ 2 Gyr and τ_c ≥ 8 Gyr, to match the selection criteria by Pfeffer et al. (2023). When looking at the bimodal fractions with respect to an age cut in the cluster populations, we find good agreement between the two models, although our simulation provides a higher number of galaxies, thus avoiding statistical fluctuations. While there appears no significant difference between star cluster populations when imposing the age constraints, there appears to be a slight difference compared to the general star cluster populations in dwarf galaxies. This result shows that the bimodal classification in dwarf galaxies can be significantly influenced by newly forming star clusters whereas the total number of star clusters in more massive galaxies far outweighs any newly formed population.

The bimodality weakly depends on the location within a galaxy. Compared to the average bimodal fraction of all star clusters we find similarity for the inner-most galaxy regions, an increase at intermediate distances, and statistical fluctuations for the outer regions, all normalised to the disk’s scale length and half-mass radius for disk- and bulge-dominated, respectively. This trend is a consequence of the location of accreted star clusters as the central regions are dominated by in situ star clusters and the outer regions by ex-situ ones. The overall trend of the bimodal fraction with galaxy stellar mass remains roughly unchanged.

Our results indicate a lower bimodality when compared to observations and the general notion that a bimodal distribution appears frequently in galaxies. To test the origin of this notion and the difference we compile a list of galaxies that were classified as uni- or bimodal from the above-mentioned literature and bin the data by galaxy mass. Overall we find bimodal fractions ranging between 50 and 100 %. However, galaxy number statistics are low. It appears likely that the discrepancy between these results and our simulations are due to differing methodologies and inhomogeneities when analysing data. For example, several studies consider a bimodality in the colour-distribution whereas others use iron abundances. Furthermore, taking the metallicity values of globular cluster candidates from Beasley et al. (2008) for NGC 5128 and applying our methodology results in a classification of uni-modal (because the weighted distance D_G is smaller than ***(73)*** $\sqrt{2}$ $Mathematical equation: $\sqrt 2 $$ ) whereas the authors classify the distribution as multi-modal. Unfortunately, the metallicities of many star clusters in galaxies are not publicly available, and we were thus not able to test the effect of different methodologies further.

Fig. 9

Fraction of galaxies exhibiting a bimodality in the metallicity distribution of their star cluster population versus stellar mass for disk-(blue) and bulge-dominated (red) galaxies. We show the bimodality when restricting the cluster population to ages τ_c ≥ 2 Gyr (light grey) and τ_c ≥ 8 Gyr (dark grey) in order to compare the results to the ones for the E-MOSAICS simulation from Pfeffer et al. (2023, dotted lines).

3.6 Caveats

Our semi-analytical approach to model star cluster populations comes with several caveats. Here we summarise a few of the most important issues.

3.6.1 Axisymmetric structures of galaxies

Several effects of non-axisymmetric features can impact star clusters, some of which are discussed below. A thorough review is presented in Renaud (2018).

Some basic assumptions of L-Galaxies potentially break down as the gas fraction of galaxies increases, resulting in a ‘clumpier’ structure (e.g. Conselice et al. 2004; Förster Schreiber et al. 2011; Shibuya et al. 2016). This raises some doubt about whether our model captures the properties of star clusters at redshifts of about one to two where clumpy galaxies start to make up a dominant fraction of the whole galaxy population (e.g. Sattari et al. 2023; Huertas-Company et al. 2024) and whether the z = 0 star cluster population evolved in a similar fashion compared to observed star clusters. However, as argued by Ono et al. (2025) the basic principle of a disk-dominated formation scenario for galaxies may be a reasonable assumption for redshifts of z ≳ 9.

A related issue is that L-Galaxies does not consider non-axisymmetric components such as bars. Molecular gas can gather at tips of the bars due to orbital crowding (Kenney & Lord 1991) leading to collisions of GMCs, triggering the formation of stars and star clusters (e.g. Davies et al. 2012; Fukui et al. 2014; Ramírez-Alegría et al. 2014). A bar can influence the dynamical evolution of star clusters as well. For instance, both Bajkova et al. (2023) and Dillamore et al. (2024) argue that the Galactic bar directly influences the orbit of a number of Milky Way globular clusters, possibly supporting accelerating their radial migration to the Galactic Centre. While halo star clusters are most likely not significantly affected by the lacking implementation of bars, objects in the galactic disk will likely be influenced, thus potentially changing the expected galactocentric distribution discussed in Section 3.4.

Finally, L-Galaxies does not consider the existence of spiral arms. As argued by Saha et al. (2010), (transient) spirals and bars can introduce additional energy sources for tidal heating that would increase the orthogonal velocity dispersion of star clusters compared to the orientation of the galactic disk. This effect could thus contribute star clusters to the halo of a galaxy without invoking galaxy-galaxy interactions, as indicated by the presence of metal-rich open clusters in the Milky Way’s halo (see Paunzen & Netopil 2006; Meibom et al. 2009; Brogaard et al. 2012; Heiter et al. 2014; Önehag et al. 2014; Straižys et al. 2014; Gustafsson et al. 2016; Hunt & Reffert 2023, for examples).

3.6.2 Impact of stellar-mass black holes

Over time, stellar mass black holes segregate towards the cluster’s centre and build up a dense core, injecting energy into the stellar system of the star cluster (e.g. Merritt et al. 2004; Mackey et al. 2008), possibly leading to the star clusters dissolution (e.g. Giersz et al. 2019). This can result in some star clusters having relatively large half-mass radii, such as Palomar 5 (Gieles et al. 2021), rivalling the most extended objects we find in our simulation. However, note that we only consider star-star interactions in the formalism of Eq. (38), thus neglecting the effect of dynamical heating due to black holes. One consequence of adding feedback from black holes is a metallicity-dependent expansion rate that was already explored in the literature (e.g. Downing 2012; Mapelli & Bressan 2013; Banerjee 2017; Chattopadhyay et al. 2022; Rostami-Shirazi et al. 2024). This may result in a change in half-mass radii of young star clusters between different galaxy morphologies that we do not detect in Figure 6. We aim to implement this feedback channel in future versions.

3.6.3 Galaxy-galaxy interactions

Interactions between galaxies result in collisions between gas clouds and can cause efficient star formation outside of galaxy disks, often including the formation of star clusters (e.g. Fellhauer & Kroupa 2005; Annibali et al. 2011; Maji et al. 2017; Randriamanakoto et al. 2019; Rodruck et al. 2023). Such star clusters, especially during the first passage of the accreted galaxy, can survive for significant time and could contribute to the z = 0 globular cluster population in the halo (e.g. Li et al. 2022). Keller et al. (2020) find that around 20% of globular clusters form during galaxy-galaxy merger events. If true, this would indicate that our model approach does not explain the origin of a significant fraction of globular clusters at z = 0. Nevertheless, as argued by the authors, repositioning of globular clusters from the dense inner-galactic regions into a galaxy’s halo is important for cluster survival, which matches our results.

4 Conclusions

We have introduced a modified version of the semi-analytical galaxy formation model ‘L-Galaxies’ (Henriques et al. 2020; Yates et al. 2021) that accounts for the formation of massive (m_c ≥ 10⁴ M_⊙) star clusters. This implementation relies on galaxy constituents to derive the bound fraction of star formation and the total star cluster mass via L-Galaxies’ prescription of star formation within galaxy disks. Star cluster masses are random realisations of an environmentally dependent CIMF that is assumed to be a truncated power-law function. Afterwards, all star clusters are assigned initial half-mass radii, metallicities, and galactocentric distances. We evolved the properties of up to 2000 individual star clusters per galaxy while taking into account the effects of stellar evolution, two-body relaxation, tidal shocks, dynamical friction, and a redistribution during galaxy mergers.

Running the simulation on output merger trees from the Millennium (Springel et al. 2005) and Millennium-II (Boylan-Kolchin et al. 2009) simulations yielded the following results.

The most massive and young (τ_c < 0.3 Gyr) star clusters in disk-dominated galaxies follow the observed empirical relationship between their absolute V-band magnitude and the total SFR of the host galaxies. There are secondary dependencies on the host galaxy’s stellar mass and cluster formation efficiency: the convolved bound fraction of star formation and the initial star cluster survival rate.
The assumption of the relationship between the sound speed of cold gas in the interstellar medium and the surface SFR directly influences the properties of young star clusters. For example, assuming that the turbulence in the interstellar medium is mainly related to gravity results in a slope value of the relationship between the V-band magnitude and the SFR being too shallow.
The star cluster mass function for young clusters exhibits a profile similar to the observational results of nearby disk-dominated galaxies (Brown & Gnedin 2021). More massive galaxies host both a larger number of star clusters and more massive ones.
The half-mass radii of star clusters are limited to about ten parsecs in a galaxy disk. In galaxy halos, depending on the galactocentric distance, star clusters can expand to half-mass radii greater than 100 pc. Our data cannot reproduce well the observational constraints of Brown & Gnedin (2021). While specific choices of the distribution of initial half-mass radii can sometimes match the observations in specific galaxy mass ranges, no prescription can match the observational constraints in all galaxy types at all galaxy masses. Reproducing these observations remains challenging for current simulations.
The mass, half-mass radius, and galactocentric distance functions of old (τ_c ≥ 6 Gyr) star clusters display complex shapes, which result from different star cluster origins (in situ, accreted, heated during galaxy mergers), and environmentally dependent prescriptions that impact star cluster evolution. Galaxy-galaxy interactions and mergers play a vital role in shaping the properties of the z = 0 star cluster population.
Our model is in excellent agreement with observations of the correlation between a bulge-dominated galaxy’s mean star cluster metallicity, as traced by iron abundance, and its host galaxy stellar mass over 4 dex. This result corroborates the importance of taking into account metal enrichment of the circum-galactic medium from supernovae, as introduced in Yates et al. (2021).
Our results predict that the presence of a bimodality in the metallicity distribution of star clusters decreases from ≈50 to ≈20% with galaxy mass at z = 0. We do not find a bimodality significantly above 50% for any galaxy morphology at any mass scale. Both different methodological approaches and the inaccessibility of a statistically significant dataset of star cluster populations do not allow for a clean comparison between our simulated galaxies and observations.

Our simulation offers a computationally efficient and flexible approach to probe different physical effects that influence the assembly history of star clusters across diverse galaxy populations in mass, type, and evolution over cosmic time. In future work we plan to look at additional aspects of the model, such as the star cluster properties of Milky Way analogues or their evolution with redshift, and to consider the formation of nuclear clusters as well as their interactions with (massive) black holes.

Data availability

The code and full data output will be published with future papers of the series. Early access may be granted upon request.

Acknowledgements

N.H. thanks F. Belfiore, I. Cabrera-Ziri, M. Gieles, J. Greene, J. Lee, M. Mapelli, D. Maschmann, T. Naab, M. Reina-Campos, A. Seth, V. Springel, and S. Torniamenti for useful discussions. N.H. was a Fellow of the International Max Planck Research School for Astronomy and Cosmic Physics at the University of Heidelberg (IMPRS-HD) and acknowledges their support. N.H. acknowledges funding through LACEGAL, a Latinamerican Chinese European Galaxy Formation Network. This project has received funding from the European Union’s HORIZON-MSCA-2021-SE-01 Research and Innovation programme under the Marie Sklodowska-Curie grant agreement number 101 086 388. N.H. acknowledges funding from the Centre national d’études spatiales (CNES) via their postdoctoral fellowship programme. S.B. acknowledges support from the Spanish Ministerio de Ciencia e Innovación through project PID2021-124243NB-C21. D.H.C. acknowledges support from the Basque Government through the Programa Predoctoral de Formación de Personal Investigador No Doctor del departamento de Ciencia, Universidades e Innovación del Gobierno Vasco. D.S. acknowledges the support by the Tsinghua Shui Mu Scholarship, the funding of the National Key R&D Program of China (grant No. 2023YFA1605600), the science research grants from the China Manned Space Project with No. CMS-CSST2021-A05, and the Tsinghua University Initiative Scientific Research Program (No. 0223080023). M.C.A. acknowledges support from Fondecyt Iniciación #11240540. M.P. acknowledges funding by the European Union (ERC, Unleash-TDEs, project number 101163093). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them.

References

Adamo, A., Östlin, G., & Zackrisson, E. 2011, MNRAS, 417, 1904 [NASA ADS] [CrossRef] [Google Scholar]
Adamo, A., Kruijssen, J. M. D., Bastian, N., Silva-Villa, E., & Ryon, J. 2015, MNRAS, 452, 246 [NASA ADS] [CrossRef] [Google Scholar]
Adamo, A., Ryon, J. E., Messa, M., et al. 2017, ApJ, 841, 131 [Google Scholar]
Adamo, A., Bradley, L. D., Vanzella, E., et al. 2024, Nature, 632, 513 [NASA ADS] [CrossRef] [Google Scholar]
Alexander, P. E. R., & Gieles, M. 2012, MNRAS, 422, 3415 [NASA ADS] [CrossRef] [Google Scholar]
Alexander, P. E. R., Gieles, M., Lamers, H. J. G. L. M., & Baumgardt, H. 2014, MNRAS, 442, 1265 [Google Scholar]
Andersson, E. P., Mac Low, M.-M., Agertz, O., Renaud, F., & Li, H. 2024, A&A, 681, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Angulo, R. E., & Hilbert, S. 2015, MNRAS, 448, 364 [CrossRef] [Google Scholar]
Angulo, R. E., & White, S. D. M. 2010, MNRAS, 405, 143 [NASA ADS] [Google Scholar]
Annibali, F., Tosi, M., Monelli, M., et al. 2009, AJ, 138, 169 [NASA ADS] [CrossRef] [Google Scholar]
Annibali, F., Tosi, M., Aloisi, A., & van der Marel, R. P. 2011, AJ, 142, 129 [NASA ADS] [CrossRef] [Google Scholar]
Ayromlou, M., Kauffmann, G., Yates, R. M., Nelson, D., & White, S. D. M. 2021, MNRAS, 505, 492 [NASA ADS] [CrossRef] [Google Scholar]
Bajkova, A. T., Smirnov, A. A., & Bobylev, V. V. 2023, AB, 78, 499 [Google Scholar]
Banerjee, S. 2017, MNRAS, 467, 524 [NASA ADS] [Google Scholar]
Banerjee, S., & Kroupa, P. 2017, A&A, 597, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Barmby, P., & Huchra, J. P. 2001, AJ, 122, 2458 [CrossRef] [Google Scholar]
Barmby, P., Huchra, J. P., Brodie, J. P., et al. 2000, AJ, 119, 727 [NASA ADS] [CrossRef] [Google Scholar]
Barrera, M., Springel, V., White, S. D. M., et al. 2023, MNRAS, 525, 6312 [NASA ADS] [CrossRef] [Google Scholar]
Bassino, L. P., & Caso, J. P. 2017, MNRAS, 466, 4259 [NASA ADS] [Google Scholar]
Bastian, N. 2008, MNRAS, 390, 759 [Google Scholar]
Bastian, N., Adamo, A., Gieles, M., et al. 2012, MNRAS, 419, 2606 [Google Scholar]
Bastian, N., Pfeffer, J., Kruijssen, J. M. D., et al. 2020, MNRAS, 498, 1050 [Google Scholar]
Baugh, C. M. 2006, Rep. Prog. Phys., 69, 3101 [NASA ADS] [CrossRef] [Google Scholar]
Baumgardt, H., & Hilker, M. 2018, MNRAS, 478, 1520 [Google Scholar]
Beasley, M. A., Bridges, T., Peng, E., et al. 2008, MNRAS, 386, 1443 [NASA ADS] [CrossRef] [Google Scholar]
Beasley, M. A., Leaman, R., Gallart, C., et al. 2019, MNRAS, 487, 1986 [Google Scholar]
Bertin, G., & Romeo, A. B. 1988, A&A, 195, 105 [NASA ADS] [Google Scholar]
Bertoldi, F., & McKee, C. F. 1992, ApJ, 395, 140 [NASA ADS] [CrossRef] [Google Scholar]
Bik, A., Lamers, H. J. G. L. M., Bastian, N., Panagia, N., & Romaniello, M. 2003, A&A, 397, 473 [CrossRef] [EDP Sciences] [Google Scholar]
Binney, J., & Tremaine, S. 2008, Galactic Dynamics, 2nd edn. (Princeton: Princeton University Press) [Google Scholar]
Blakeslee, J. P. 1999, AJ, 118, 1506 [NASA ADS] [CrossRef] [Google Scholar]
Blom, C., Spitler, L. R., & Forbes, D. A. 2012, MNRAS, 420, 37 [Google Scholar]
Bluck, A. F. L., Mendel, T. J., Ellison, S. L., et al. 2016, MNRAS, 462, 2559 [Google Scholar]
Bonoli, S., Marulli, F., Springel, V., et al. 2009, MNRAS, 396, 423 [NASA ADS] [CrossRef] [Google Scholar]
Bonoli, S., Mayer, L., & Callegari, S. 2014, MNRAS, 437, 1576 [NASA ADS] [CrossRef] [Google Scholar]
Bottema, R. 1993, A&A, 275, 16 [NASA ADS] [Google Scholar]
Boylan-Kolchin, M., Springel, V., White, S. D. M., Jenkins, A., & Lemson, G. 2009, MNRAS, 398, 1150 [Google Scholar]
Bradley, L. D., Adamo, A., Vanzella, E., et al. 2024, ApJ, 22 [Google Scholar]
Brodie, J. P., & Strader, J. 2006, ARAA, 44, 193 [Google Scholar]
Brodie, J. P., Usher, C., Conroy, C., et al. 2012, ApJ, 759, L33 [NASA ADS] [CrossRef] [Google Scholar]
Brogaard, K., VandenBerg, D. A., Bruntt, H., et al. 2012, A&A, 543, A106 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Brown, G., & Gnedin, O. Y. 2021, MNRAS, 508, 5935 [NASA ADS] [CrossRef] [Google Scholar]
Brown, G., & Gnedin, O. Y. 2022, MNRAS, 514, 280 [Google Scholar]
Burkert, A., & Forbes, D. A. 2020, AJ, 159, 56 [Google Scholar]
Burkhart, B. 2018, ApJ, 863, 118 [CrossRef] [Google Scholar]
Caldwell, N., & Romanowsky, A. J. 2016, ApJ, 824, 42 [NASA ADS] [CrossRef] [Google Scholar]
Callingham, T. M., Cautun, M., Deason, A. J., et al. 2022, MNRAS, 513, 4107 [NASA ADS] [CrossRef] [Google Scholar]
Calura, F., Pascale, R., Agertz, O., et al. 2025, A&A, 698, A207 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Calzetti, D., Lee, J. C., Sabbi, E., et al. 2015, AJ, 149, 51 [Google Scholar]
Cantiello, M., Capaccioli, M., Napolitano, N., et al. 2015, A&A, 576, A14 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Carretta, E., Gratton, R. G., Clementini, G., & Fusi Pecci, F. 2000, ApJ, 533, 215 [NASA ADS] [CrossRef] [Google Scholar]
Chabrier, G. 2003, PASP, 115, 763 [Google Scholar]
Chandrasekhar, S. 1943, ApJ, 97, 255 [Google Scholar]
Chattopadhyay, D., Hurley, J., Stevenson, S., & Raidani, A. 2022, MNRAS, 513, 4527 [NASA ADS] [CrossRef] [Google Scholar]
Chen, Y., Mo, H., & Wang, H. 2025, MNRAS, 540, 1235 [Google Scholar]
Chevance, M., Kruijssen, J. M. D., Hygate, A. P. S., et al. 2020, MNRAS, 493, 2872 [NASA ADS] [CrossRef] [Google Scholar]
Choksi, N., & Gnedin, O. Y. 2019, MNRAS, 486, 331 [NASA ADS] [CrossRef] [Google Scholar]
Choksi, N., & Kruijssen, J. M. D. 2021, MNRAS, 507, 5492 [Google Scholar]
Choksi, N., Behroozi, P., Volonteri, M., et al. 2017, MNRAS, 472, 1526 [Google Scholar]
Claeyssens, A., Adamo, A., Richard, J., et al. 2023, MNRAS, 520, 2180 [NASA ADS] [CrossRef] [Google Scholar]
Cohen, J. G., Blakeslee, J. P., & Ryzhov, A. 1998, ApJ, 496, 808 [NASA ADS] [CrossRef] [Google Scholar]
Colombo, D., Hughes, A., Schinnerer, E., et al. 2014, ApJ, 784, 3 [NASA ADS] [CrossRef] [Google Scholar]
Conroy, C., & Gunn, J. E. 2010a, Astrophysics Source Code Library [record ascl:1010.043] [Google Scholar]
Conroy, C., & Gunn, J. E. 2010b, ApJ, 712, 833 [Google Scholar]
Conroy, C., White, M., & Gunn, J. E. 2010, ApJ, 708, 58 [NASA ADS] [CrossRef] [Google Scholar]
Conroy, C., Gunn, J. E., & White, M. 2009, ApJ, 699, 486 [Google Scholar]
Conselice, C. J., Grogin, N. A., Jogee, S., et al. 2004, ApJ, 600, L139 [NASA ADS] [CrossRef] [Google Scholar]
Cook, D. O., Seth, A. C., Dale, D. A., et al. 2012, ApJ, 751, 100 [NASA ADS] [CrossRef] [Google Scholar]
Cook, D. O., Lee, J. C., Adamo, A., et al. 2023, MNRAS, 519, 3749 [CrossRef] [Google Scholar]
Correnti, M., Gennaro, M., Kalirai, J. S., Brown, T. M., & Calamida, A. 2016, ApJ, 823, 18 [Google Scholar]
Crain, R. A., Schaye, J., Bower, R. G., et al. 2015, MNRAS, 450, 1937 [NASA ADS] [CrossRef] [Google Scholar]
Croton, D. J. 2006, MNRAS, 369, 1808 [NASA ADS] [CrossRef] [Google Scholar]
Croton, D. J., Springel, V., White, S. D. M., et al. 2006, MNRAS, 365, 11 [Google Scholar]
da Silva, R. L., Fumagalli, M., & Krumholz, M. 2012, ApJ, 745, 145 [CrossRef] [Google Scholar]
da Silva, R. L., Fumagalli, M., & Krumholz, M. 2014, MNRAS, 444, 3275 [Google Scholar]
Davies, B., de La Fuente, D., Najarro, F., et al. 2012, MNRAS, 419, 1860 [Google Scholar]
Davis, M., Efstathiou, G., Frenk, C. S., & White, S. D. M. 1985, ApJ, 292, 371 [Google Scholar]
De Lucia, G., & Blaizot, J. 2007, MNRAS, 375, 2 [Google Scholar]
De Lucia, G., Tornatore, L., Frenk, C. S., et al. 2014, MNRAS, 445, 970 [Google Scholar]
De Lucia, G., Kruijssen, J. M. D., Trujillo-Gomez, S., Hirschmann, M., & Xie, L. 2024, MNRAS, 530, 2760 [NASA ADS] [CrossRef] [Google Scholar]
del P. Lagos, C., Bravo, M., Tobar, R., et al. 2024, MNRAS, 531, 3551 [Google Scholar]
Dillamore, A. M., Monty, S., Belokurov, V., & Evans, N. W. 2024, ApJ, 971, L4 [Google Scholar]
Dobbs, C. L., Burkert, A., & Pringle, J. E. 2011, MNRAS, 417, 1318 [NASA ADS] [CrossRef] [Google Scholar]
Dolag, K., Borgani, S., Murante, G., & Springel, V. 2009, MNRAS, 399, 497 [Google Scholar]
Dolfi, A., Pfeffer, J., Forbes, D. A., et al. 2022, MNRAS, 511, 3179 [Google Scholar]
Downing, J. M. B. 2012, MNRAS, 425, 2234 [Google Scholar]
Dutton, A. A., & Macció, A. V. 2014, MNRAS, 441, 3359 [NASA ADS] [CrossRef] [Google Scholar]
Ekström, S., Georgy, C., Eggenberger, P., et al. 2012, A&A, 537, A146 [Google Scholar]
El-Badry, K., Quataert, E., Weisz, D. R., Choksi, N., & Boylan-Kolchin, M. 2019, MNRAS, 482, 4528 [Google Scholar]
Elmegreen, B., & Lahén, N. 2024, OJA, 7, 68 [Google Scholar]
Elmegreen, D. M., Elmegreen, B. G., Ravindranath, S., & Coe, D. A. 2007, ApJ, 658, 763 [Google Scholar]
Emig, K. L., Bolatto, A. D., Leroy, A. K., et al. 2020, ApJ, 903, 50 [NASA ADS] [CrossRef] [Google Scholar]
Emsellem, E., Schinnerer, E., Santoro, F., et al. 2022, A&A, 659, A191 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Escala, A., Larson, R. B., Coppi, P. S., & Mardones, D. 2004, ApJ, 607, 765 [NASA ADS] [CrossRef] [Google Scholar]
Escudero, C. G., Faifer, F. R., Bassino, L. P., Calderón, J. P., & Caso, J. P. 2015, MNRAS, 449, 612 [Google Scholar]
Fahrion, K., Lyubenova, M., Hilker, M., et al. 2020, A&A, 637, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Fahrion, K., Lyubenova, M., van de Ven, G., et al. 2021, A&A, 650, A137 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Fellhauer, M., & Kroupa, P. 2005, MNRAS, 359, 223 [Google Scholar]
Feng, C.-C., Lin, L.-H., Wang, H.-H., & Taam, R. E. 2014, ApJ, 785, 103 [Google Scholar]
Filion, C., McClure, R. L., Weinberg, M. D., D’Onghia, E., & Daniel, K. J. 2023, MNRAS, 524, 276 [Google Scholar]
Flynn, C., Holmberg, J., Portinari, L., Fuchs, B., & Jahreiß, H. 2006, MNRAS, 372, 1149 [NASA ADS] [CrossRef] [Google Scholar]
Forbes, D. A., & Bridges, T. 2010, MNRAS, 404, 1203 [NASA ADS] [Google Scholar]
Forbes, D. A., & Romanowsky, A. J. 2023, MNRAS, 520, L58 [Google Scholar]
Forbes, D. A., Bastian, N., Gieles, M., et al. 2018, Proc. R. Soc. London Ser. A, 474, 20170616 [Google Scholar]
Förster Schreiber, N. M., Shapley, A. E., Erb, D. K., et al. 2011, ApJ, 731, 65 [CrossRef] [Google Scholar]
Freeman, K. C. 1970, ApJ, 160, 811 [Google Scholar]
Fu, J., Guo, Q., Kauffmann, G., & Krumholz, M. R. 2010, MNRAS, 409, 515 [Google Scholar]
Fu, J., Kauffmann, G., Li, C., & Qi, G. 2012, MNRAS, 424, 2701 [Google Scholar]
Fu, J., Kauffmann, G., Huang, M.-L., et al. 2013, MNRAS, 434, 1531 [NASA ADS] [CrossRef] [Google Scholar]
Fukui, Y., Ohama, A., Hanaoka, N., et al. 2014, ApJ, 780, 36 [Google Scholar]
Gabrielpillai, A., Somerville, R. S., Genel, S., et al. 2022, MNRAS, 517, 6091 [Google Scholar]
Garcia, F. A. B., Ricotti, M., Sugimura, K., & Park, J. 2023, MNRAS, 522, 2495 [NASA ADS] [CrossRef] [Google Scholar]
Garro, E. R., Minniti, D., & Fernández-Trincado, J. G. 2024, A&A, 687, A214 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Genzel, R., Newman, S., Jones, T., et al. 2011, ApJ, 733, 101 [Google Scholar]
Georgy, C., Ekström, S., Eggenberger, P., et al. 2013, A&A, 558, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Gieles, M., & Renaud, F. 2016, MNRAS, 463, L103 [Google Scholar]
Gieles, M., Heggie, D. C., & Zhao, H. 2011, MNRAS, 413, 2509 [Google Scholar]
Gieles, M., Erkal, D., Antonini, F., Balbinot, E., & Peñarrubia, J. 2021, NatAs, 5, 957 [Google Scholar]
Giersz, M., & Heggie, D. C. 1994, MNRAS, 268, 257 [NASA ADS] [CrossRef] [Google Scholar]
Giersz, M., Askar, A., Wang, L., et al. 2019, MNRAS, 487, 2412 [NASA ADS] [CrossRef] [Google Scholar]
Goddard, Q. E., Bastian, N., & Kennicutt, R. C. 2010, MNRAS, 405, 857 [NASA ADS] [Google Scholar]
Green, A. W., Glazebrook, K., McGregor, P. J., et al. 2014, MNRAS, 437, 1070 [Google Scholar]
Guo, S., & Qi, Z. 2023, Universe, 9, 252 [Google Scholar]
Guo, Q., White, S., Boylan-Kolchin, M., et al. 2011, MNRAS, 413, 101 [Google Scholar]
Gustafsson, B., Church, R. P., Davies, M. B., & Rickman, H. 2016, A&A, 593, A85 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Hackshaw, Z., Hawkins, K., Filion, C., et al. 2024, ApJ, 255, 19 [Google Scholar]
Hammer, F., Li, H., Mamon, G. A., et al. 2023, MNRAS, 519, 5059 [CrossRef] [Google Scholar]
Hammer, F., Wang, J., Mamon, G. A., et al. 2024, MNRAS, 527, 2718 [Google Scholar]
Harris, W. E., Ciccone, S. M., Eadie, G. M., et al. 2017, ApJ, 835, 101 [NASA ADS] [CrossRef] [Google Scholar]
Hawkins, K. 2023, MNRAS, 525, 3318 [NASA ADS] [CrossRef] [Google Scholar]
Heiter, U., Soubiran, C., Netopil, M., & Paunzen, E. 2014, A&A, 561, A93 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Hénon, M. 1961, Annales d’Astrophysique, 24, 369 [NASA ADS] [Google Scholar]
Henriques, B. M. B., White, S. D. M., Thomas, P. A., et al. 2015, MNRAS, 451, 2663 [Google Scholar]
Henriques, B. M. B., Yates, R. M., Fu, J., et al. 2020, MNRAS, 491, 5795 [NASA ADS] [CrossRef] [Google Scholar]
Heyer, M. H., Corbelli, E., Schneider, S. E., & Young, J. S. 2004, ApJ, 602, 723 [NASA ADS] [CrossRef] [Google Scholar]
Heyer, M., Krawczyk, C., Duval, J., & Jackson, J. M. 2009, ApJ, 699, 1092 [Google Scholar]
Hirschmann, M., De Lucia, G., & Fontanot, F. 2016, MNRAS, 461, 1760 [Google Scholar]
Hislop, J. M., Naab, T., Steinwandel, U. P., et al. 2022, MNRAS, 509, 5938 [Google Scholar]
Hixenbaugh, K., Chandar, P., & Mok, A. 2022, AJ, 163, 271 [Google Scholar]
Ho, I.-T., Seibert, M., Meidt, S. E., et al. 2017, ApJ, 846, 39 [NASA ADS] [CrossRef] [Google Scholar]
Ho, I.-T., Kreckel, K., Meidt, S. E., et al. 2019, ApJ, 885, L31 [NASA ADS] [CrossRef] [Google Scholar]
Hollyhead, K., Adamo, A., Bastian, N., Gieles, M., & Ryon, J. E. 2016, MNRAS, 460, 2087 [Google Scholar]
Hopkins, P. F., Quataert, E., & Murray, N. 2012, MNRAS, 421, 3488 [NASA ADS] [CrossRef] [Google Scholar]
Horta, D., Hughes, M. E., Pfeffer, J. L., et al. 2021, MNRAS, 500, 4768 [Google Scholar]
Hoyer, N., Neumayer, N., Georgiev, I. Y., Seth, A. C., & Greene, J. E. 2021, MNRAS, 507, 3246 [NASA ADS] [CrossRef] [Google Scholar]
Huchra, J. P., Brodie, J. P., & Kent, S. M. 1991, ApJ, 370, 495 [Google Scholar]
Hudson, M. J., Harris, G. L., & Harris, W. E. 2014, ApJ, 787, L5 [Google Scholar]
Huertas-Company, M., Iyer, K. G., Angeloudi, E., et al. 2024, A&A, 685, A48 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Hughes, M. E., Pfeffer, J., Martig, M., et al. 2019, MNRAS, 482, 2795 [Google Scholar]
Hughes, M. E., Pfeffer, J. L., Martig, M., et al. 2020, MNRAS, 491, 4012 [NASA ADS] [CrossRef] [Google Scholar]
Hughes, M. E., Pfeffer, J. L., Bastian, N., et al. 2022, MNRAS, 510, 6190 [Google Scholar]
Hunt, E. L., & Reffert, S. 2023, A&A, 673, A114 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Hunter, D. A., Elmegreen, B. G., Dupuy, T. J., & Mortonson, M. 2003, AJ, 126, 1836 [NASA ADS] [CrossRef] [Google Scholar]
Huxor, A. P., Mackey, A. D., Ferguson, A. M. N., et al. 2014, MNRAS, 442, 2165 [Google Scholar]
Hwang, H.-C., Barrera-Ballesteros, J. K., Heckman, T. M., et al. 2019, ApJ, 872, 144 [NASA ADS] [CrossRef] [Google Scholar]
Ines Ennis, A., Caso, J. P., & Bassino, L. P. 2024, OJA, 12 [Google Scholar]
Irodotou, D., Thomas, P. A., Henriques, B. M. B., Sargent, M. T., & Hislop, J. M. 2019, MNRAS, 489, 3609 [Google Scholar]
Izquierdo-Villalba, D., Bonoli, S., Spinoso, D., et al. 2019, MNRAS, 488, 609 [NASA ADS] [CrossRef] [Google Scholar]
Izquierdo-Villalba, D., Bonoli, S., Dotti, M., et al. 2020, MNRAS, 495, 4681 [NASA ADS] [CrossRef] [Google Scholar]
Izquierdo-Villalba, D., Sesana, A., Bonoli, S., & Colpi, M. 2022, MNRAS, 509, 3488 [Google Scholar]
Jaffe, W. 1983, MNRAS, 202, 995 [NASA ADS] [Google Scholar]
Johnson, B., Foreman-Mackey, D., Sick, J., et al. 2023, https://doi.org/10.5281/zenodo.10026684 [Google Scholar]
Johnson, K. E., Leitherer, C., Vacca, W. D., & Conti, P. S. 2000, AJ, 120, 1273 [NASA ADS] [CrossRef] [Google Scholar]
Kacharov, N., Neumayer, N., Seth, A. C., et al. 2018, MNRAS, 480, 1973 [Google Scholar]
Karam, J., & Sills, A. 2022, MNRAS, 513, 6095 [Google Scholar]
Karam, J., & Sills, A. 2023, MNRAS, 521, 5557 [Google Scholar]
Katz, H., & Ricotti, M. 2013, MNRAS, 432, 3250 [Google Scholar]
Kauffmann, G., & Haehnelt, M. 2000, MNRAS, 311, 576 [Google Scholar]
Kaviraj, S., Ferreras, I., Yoon, S.-J., & Yi, S. K. 2005, A&A, 439, 913 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Keller, B. W., Kruijssen, J. M. D., Pfeffer, J., et al. 2020, MNRAS, 495, 4248 [NASA ADS] [CrossRef] [Google Scholar]
Kenney, J. D. P., & Lord, S. D. 1991, ApJ, 381, 118 [NASA ADS] [CrossRef] [Google Scholar]
Kennicutt, R. C. 1998, ARAA, 36, 189 [NASA ADS] [CrossRef] [Google Scholar]
Kewley, L. J., & Ellison, S. L. 2008, ApJ, 681, 1183 [Google Scholar]
Kikuchihara, S., Ouchi, M., Ono, Y., et al. 2020, ApJ, 893, 60 [NASA ADS] [CrossRef] [Google Scholar]
King, I. 1962, AJ, 67, 471 [Google Scholar]
Krauss, L. M., & Chaboyer, B. 2003, Science, 229, 65 [Google Scholar]
Kravtsov, A. V., & Gnedin, O. Y. 2005, ApJ, 623, 650 [Google Scholar]
Kreckel, K., Ho, I.-T., Blanc, G. A., et al. 2019, ApJ, 887, 80 [NASA ADS] [CrossRef] [Google Scholar]
Kritsuk, A. G., Norman, M. L., Padoan, P., & Wagner, R. 2007, ApJ, 665, 416 [NASA ADS] [CrossRef] [Google Scholar]
Kritsuk, A. G., Ustyugov, S. D., & Norman, M. L. 2017, New J. Phys., 19, 065003 [Google Scholar]
Kruijssen, J. M. D. 2012, MNRAS, 426, 3008 [Google Scholar]
Kruijssen, J. M. D. 2014, Class. Quantum Gravity, 31, 244006 [NASA ADS] [CrossRef] [Google Scholar]
Kruijssen, J. M. D. 2015, MNRAS, 454, 1658 [Google Scholar]
Kruijssen, J. M. D., & Cooper, A. P. 2012, MNRAS, 420, 340 [Google Scholar]
Kruijssen, J. M. D., Pfeffer, J. L., Crain, R. A., & Bastian, N. 2019a, MNRAS, 486, 3134 [NASA ADS] [CrossRef] [Google Scholar]
Kruijssen, J. M. D., Pfeffer, J. L., Reina-Campos, M., Crain, R. A., & Bastian, N. 2019b, MNRAS, 486, 3180 [Google Scholar]
Kruijssen, J. M. D., Pfeffer, J. L., Chevance, M., et al. 2020, MNRAS, 498, 2472 [NASA ADS] [CrossRef] [Google Scholar]
Krumholz, M. R., & Burkhart, B. 2016, MNRAS, 458, 1671 [NASA ADS] [CrossRef] [Google Scholar]
Krumholz, M. R., & McKee, C. F. 2005, ApJ, 630, 260 [Google Scholar]
Krumholz, M. R., McKee, C. F., & Tumlinson, J. 2009, ApJ, 693, 216 [NASA ADS] [CrossRef] [Google Scholar]
Krumholz, M. R., Fumagalli, M., da Silva, R. L., Rendahl, T., & Parra, J. 2015, MNRAS, 452, 1447 [NASA ADS] [CrossRef] [Google Scholar]
Krumholz, M. R., Burkhart, B., Forbes, J. C., & Crocker, R. M. 2018, MNRAS, 477, 2716 [CrossRef] [Google Scholar]
Kuijken, K., & Gilmore, G. 1989, MNRAS, 239, 571 [NASA ADS] [CrossRef] [Google Scholar]
Kundu, A., & Whitmore, B. C. 2001, AJ, 121, 2950 [NASA ADS] [CrossRef] [Google Scholar]
Lahén, N., Naab, T., Johansson, P. H., et al. 2020, ApJ, 891, 2 [CrossRef] [Google Scholar]
Lahén, N., Naab, T., Kauffmann, G., et al. 2023, MNRAS, 522, 3092 [CrossRef] [Google Scholar]
Lahén, N., Naab, T., & Szécsi, D. 2024, MNRAS, 530, 645 [CrossRef] [Google Scholar]
Lahén, N., Rantala, A., Naab, T., et al. 2025, MNRAS, 538, 2129 [Google Scholar]
Lake, W., Noaz, S., Chiou, Y. S., et al. 2021, ApJ, 922, 86 [Google Scholar]
Lake, W., Naoz, S., Marinacci, F., et al. 2023, ApJ, 956, L7 [Google Scholar]
Lamers, H. J. G. L. M., Gieles, M., & Portegies Zwart, S. F. 2005, A&A, 429, 173 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Larsen, S. S. 2002, AJ, 124, 1393 [NASA ADS] [CrossRef] [Google Scholar]
Larsen, S. S. 2010, Philos. Trans. R. Soc. A, 368, 867 [Google Scholar]
Larsen, S. S., Romanowsky, A. J., Brodie, J. P., & Wasserman, A. 2020, Science, 370, 970 [NASA ADS] [CrossRef] [Google Scholar]
Larson, R. B. 1981, MNRAS, 194, 809 [Google Scholar]
Le, M. N., & Cooper, A. P. 2025, ApJ, 978, 33 [Google Scholar]
Leaman, R., VandenBerg, D. A., & Mendel, J. T. 2013, MNRAS, 436, 122 [Google Scholar]
Lee, J. C., Sandstrom, K. M., Leroy, A. K., et al. 2023, ApJ, 944, L17 [NASA ADS] [CrossRef] [Google Scholar]
Lee, J. C., Whitmore, B. C., Thilker, D. A., et al. 2022, ApJSS, 258, 10 [Google Scholar]
Lee, M. G. 2003, JKAS, 36, 189 [Google Scholar]
Lehnert, M. D., Nesvadba, N. P. H., Le Tiran, L., et al. 2009, ApJ, 699, 1660 [NASA ADS] [CrossRef] [Google Scholar]
Leroy, A. K., Walter, F., Brinks, E., et al. 2008, AJ, 136, 2782 [Google Scholar]
Levy, R. C., Bolatto, A. D., Mayya, D., et al. 2024, ApJ, 973, 16 [Google Scholar]
Li, H., & Gnedin, O. Y. 2019, MNRAS, 486, 4030 [NASA ADS] [CrossRef] [Google Scholar]
Li, H., Gnedin, O. Y., & Gnedin, N. Y. 2018, ApJ, 861, 107 [Google Scholar]
Li, H., Gnedin, O. Y., Gnedin, N. Y., et al. 2017, ApJ, 834, 69 [NASA ADS] [CrossRef] [Google Scholar]
Li, H., Vogelberger, M., Bryan, G. L., et al. 2022, MNRAS, 514, 265 [Google Scholar]
Lim, S., & Lee, M. G. 2015, ApJ, 804, 123 [NASA ADS] [CrossRef] [Google Scholar]
Lin, C. C., & Shu, F. H. 1966, PNAS, 55, 229 [CrossRef] [Google Scholar]
Lomelí-Núñez, L., Mayya, Y. D., Rodríguez-Merino, L. H., et al. 2024, MNRAS, 528, 1445 [Google Scholar]
Mackey, A. D., Wilkinson, M. I., Davies, M. B., & Gilmore, G. F. 2008, MNRAS, 386, 65 [NASA ADS] [CrossRef] [Google Scholar]
Mai, Y., Croom, S. M., Wisnioski, E., et al. 2024, MNRAS, 533, 15 [Google Scholar]
Maji, M., Zhu, Q., Li, Y., et al. 2017, ApJ, 844, 108 [NASA ADS] [CrossRef] [Google Scholar]
Mapelli, M., & Bressan, A. 2013, MNRAS, 430, 3120 [Google Scholar]
Maraston, C., Bastian, N., Saglia, R. P., et al. 2004, A&A, 416, 467 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Martin, N. F., Venn, K. A., Aguado, D. S., et al. 2022, Nature, 601, 45 [NASA ADS] [CrossRef] [Google Scholar]
Maschmann, D., Lee, J. C., Thilker, D. A., et al. 2024, ApJSS, 273, 14 [Google Scholar]
Matzner, C. D., & McKee, C. F. 2000, ApJ, 545, 364 [NASA ADS] [CrossRef] [Google Scholar]
McCrady, N., & Graham, J. R. 2007, ApJ, 663, 844 [NASA ADS] [Google Scholar]
McKee, C. F., & Krumholz, M. R. 2010, ApJ, 709, 308 [NASA ADS] [CrossRef] [Google Scholar]
McKee, C. F., & Tan, J. C. 2003, ApJ, 585, 850 [Google Scholar]
Meibom, S., Grundahl, F., Clausen, J. V., et al. 2009, AJ, 137, 5086 [Google Scholar]
Merritt, D., Piatek, S., Portegies Zwart, S., & Hemsendorf, M. 2004, ApJ, 608, L25 [NASA ADS] [CrossRef] [Google Scholar]
Metha, B., Trenti, M., & Chu, T. 2021, MNRAS, 508, 489 [NASA ADS] [CrossRef] [Google Scholar]
Metha, B., Trenti, M., Battisti, A., & Chu, T. 2024, MNRAS, 529, 104 [Google Scholar]
Meynet, G., Maeder, A., Schaller, G., Schaerer, D., & Charbonnel, C. 1994, A&ASS, 103, 97 [Google Scholar]
Minniti, D., Geisler, D., Alonso-García, J., et al. 2017, ApJ, 849, L24 [CrossRef] [Google Scholar]
Monaco, P., Fontanot, F., & Taffoni, G. 2007, MNRAS, 375, 1189 [NASA ADS] [CrossRef] [Google Scholar]
Mowla, L., Iyer, K. G., Desprez, G., et al. 2022, ApJ, 937, L35 [NASA ADS] [CrossRef] [Google Scholar]
Mowla, L., Iyer, K., Asada, Y., et al. 2024, Nature, 636, 332 [Google Scholar]
Muratov, A. L., & Gnedin, O. Y. 2010, ApJ, 718, 1266 [Google Scholar]
Murphy, G. G., Yates, R. M., & Mohamed, S. S. 2022, MNRAS, 510, 1945 [Google Scholar]
Myeong, G. C., Evans, N. W., Belokurov, V., Sanders, L. J., & Koposov, S. E. 2018, ApJ, 863, L28 [NASA ADS] [CrossRef] [Google Scholar]
Myers, P. C., & Goodman, A. A. 1988, ApJ, 329, 392 [Google Scholar]
Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563 [Google Scholar]
Newton, O., Davies, J. J., Pfeffer, J., et al. 2025, MNRAS, 542, 591 [Google Scholar]
Nordlund, Å., & Padoan, P. 1999, in Interstellar Turbulence, eds. J. Franco, & A. Carraminana (Cambridge: Cambridge University Press), 218 [Google Scholar]
Ohlson, D., Seth, A. C., Gallo, E., Baldassare, V. F., & Greene, J. E. 2023, AJ, 23 [Google Scholar]
Oklopčić, A., Hopkins, P. F., Feldmann, R., et al. 2017, MNRAS, 465, 952 [CrossRef] [Google Scholar]
Önehag, A., Gustafsson, B., & Korn, A. 2014, A&A, 562, A102 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Ono, Y., Ouchi, M., Harikane, Y., et al. 2025, ApJ, 991, 222 [Google Scholar]
Ostriker, E. C. 1999, ApJ, 513, 252 [Google Scholar]
Ostriker, E. C., Stone, J. M., & Gammie, C. F. 2001, ApJ, 546, 980 [Google Scholar]
Pace, A. B. 2025, OJA, 8, 142 [Google Scholar]
Padoan, P., & Nordlund, Å. 2002, ApJ, 576, 870 [NASA ADS] [CrossRef] [Google Scholar]
Padoan, P., & Nordlund, Å. 2011, ApJ, 730, 40 [Google Scholar]
Padoan, P., Jimenez, R., & Jones, B. 1997, MNRAS, 185, 711 [Google Scholar]
Parente, M., Ragone-Figueroa, C., Granato, G. L., & Lapi, A. 2023, MNRAS, 521, 6105 [NASA ADS] [CrossRef] [Google Scholar]
Parente, M., Ragone-Figueroa, C., López, P., et al. 2024, ApJ, 966, 154 [NASA ADS] [CrossRef] [Google Scholar]
Pasquali, A., Bik, A., Zibetti, S., et al. 2011, AJ, 141, 132 [NASA ADS] [CrossRef] [Google Scholar]
Pastorello, N., Forbes, D. A., Usher, C., et al. 2015, MNRAS, 451, 2625 [Google Scholar]
Paunzen, E., & Netopil, M. 2006, MNRAS, 371, 1641 [NASA ADS] [CrossRef] [Google Scholar]
Pechetti, R., Seth, A. C., Neumayer, N., et al. 2020, ApJ, 900, 32 [NASA ADS] [CrossRef] [Google Scholar]
Peng, E. W., Jordán, A., Côté, P., et al. 2006, ApJ, 639, 95 [Google Scholar]
Peng, E. W., Jordán, A., Côté, P., et al. 2008, ApJ, 681, 197 [NASA ADS] [CrossRef] [Google Scholar]
Perrett, K. M., Bridges, T. J., Hanes, D. A., et al. 2002, AJ, 123, 2490 [Google Scholar]
Pfeffer, J., Kruijssen, J. M. D., Crain, R. A., & Bastian, N. 2018, MNRAS, 475, 4309 [Google Scholar]
Pfeffer, J., Bastian, N., Crain, R. A., et al. 2019a, MNRAS, 487, 4550 [Google Scholar]
Pfeffer, J., Bastian, N., Kruijssen, J. M. D., et al. 2019b, MNRAS, 490, 1714 [Google Scholar]
Pfeffer, J., Kruijssen, J. M. D., Bastian, N., Crain, R. A., & Trujillo-Gomez, S. 2023, MNRAS, 519, 5384 [Google Scholar]
Pfeffer, J., Janssens, S. R., Buzzo, M. L., et al. 2024, MNRAS, 429, 4914 [Google Scholar]
Pfeffer, J., Forbes, D. A., Romanowsky, A. J., et al. 2025, MNRAS, 536, 1878 [Google Scholar]
Planck Collaboration XVI. 2014, A&A, 571, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Plummer, H. C. 1911, MNRAS, 71, 460 [Google Scholar]
Poggio, E., Recio-Blano, A., Palicio, P. A., et al. 2022, A&A, 666, L4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Polkas, M., Bonoli, S., Bortolas, E., et al. 2024, A&A, 689, A204 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Portegies Zwart, S. F., McMillan, S. L. W., & Gieles, M. 2010, ARAA, 48, 431 [Google Scholar]
Prieto, J. L., & Gnedin, O. Y. 2008, ApJ, 689, 919 [Google Scholar]
Rafelski, M., & Zaritsky, D. 2005, AJ, 129, 2701 [NASA ADS] [CrossRef] [Google Scholar]
Rafikov, R. R. 2001, MNRAS, 323, 445 [NASA ADS] [CrossRef] [Google Scholar]
Ramírez-Alegría, S., Borissova, J., Chené, A. B., et al. 2014, A&A, 564, L9 [Google Scholar]
Randriamanakoto, Z., Väisänen, P., Ryder, S. D., & Ranaivomanana, P. 2019, MNRAS, 482, 2530 [Google Scholar]
Reina-Campos, M., & Kruijssen, J. M. D. 2017, MNRAS, 469, 1282 [CrossRef] [Google Scholar]
Reina-Campos, M., Kruijssen, J. M. D., Pfeffer, J. L., Bastian, N., & Crain, R. A. 2019, MNRAS, 486, 5838 [NASA ADS] [CrossRef] [Google Scholar]
Reina-Campos, M., Hughes, M. E., Kruijssen, J. M. D., et al. 2020, MNRAS, 493, 3422 [Google Scholar]
Reina-Campos, M., Keller, B. W., Kruijssen, J. M. D., et al. 2022a, MNRAS, 517, 3144 [NASA ADS] [CrossRef] [Google Scholar]
Reina-Campos, M., Trujillo-Gomez, S., Deason, A. J., et al. 2022b, MNRAS, 513, 3925 [NASA ADS] [CrossRef] [Google Scholar]
Reina-Campos, M., Sills, A., & Bichon, G. 2023a, MNRAS, 524, 968 [NASA ADS] [CrossRef] [Google Scholar]
Reina-Campos, M., Trujillo-Gomez, S., Pfeffer, J. L., et al. 2023b, MNRAS, 521, 6368 [NASA ADS] [CrossRef] [Google Scholar]
Reina-Campos, M., Gnedin, O. Y., Sills, A., & Li, H. 2025, ApJ, 978, 15 [Google Scholar]
Renaud, F. 2018, New Astron. Rev., 81, 1 [CrossRef] [Google Scholar]
Renaud, F. 2020, IAU Symp., 351, 40 [Google Scholar]
Renaud, F., Gieles, M., & Boily, C. M. 2011, MNRAS, 418, 759 [NASA ADS] [CrossRef] [Google Scholar]
Rodruck, M., Charlton, J., Borthakur, S., et al. 2023, MNRAS, 526, 2341 [Google Scholar]
Romeo, A. B., & Wiegert, J. 2011, MNRAS, 416, 1191 [NASA ADS] [CrossRef] [Google Scholar]
Rosolowsky, E., & Blitz, L. 2005, ApJ, 623, 826 [NASA ADS] [CrossRef] [Google Scholar]
Rostami-Shirazi, A., Zonoozi, A. H., Haghi, H., & Rabiee, M. 2024, MNRAS, 535, 3489 [Google Scholar]
Ryon, J. E., Adamo, A., Bastian, N., et al. 2014, AJ, 148, 33 [Google Scholar]
Ryon, J. E., Bastian, N., Adamo, A., et al. 2015, MNRAS, 452, 525 [CrossRef] [Google Scholar]
Ryon, J. E., Gallagher, J. S., Smith, L. J., et al. 2017, ApJ, 841, 92 [Google Scholar]
Safronov, V. S. 1960, Annales d’Astrophysique, 23, 979 [NASA ADS] [Google Scholar]
Saha, K., tseng, Y.-H., & Taam, R. E. 2010, ApJ, 721, 1878 [Google Scholar]
Salmon, B., Coe, D., Bradley, L., et al. 2018, ApJ, 864, L22 [Google Scholar]
Sánchez-Menguiano, L., Sánchez, S. F., Kawata, D., et al. 2016, ApJ, 830, L40 [CrossRef] [Google Scholar]
Sánchez-Menguiano, L., Sánchez, S. F., Pérez, I., et al. 2018, A&A, 609, A119 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Sanders, R. L., Shapley, A. E., Jones, T., et al. 2021, ApJ, 914, 19 [CrossRef] [Google Scholar]
Sattari, Z., Mobasher, B., Chartab, N., et al. 2023, ApJ, 951, 147 [NASA ADS] [CrossRef] [Google Scholar]
Schaller, G., Schaerer, D., Meynet, G., & Maeder, A. 1992, A&ASS, 96, 269 [Google Scholar]
Schaye, J., Crain, R. A., Bower, R. G., et al. 2015, MNRAS, 446, 521 [Google Scholar]
Schuster, K. F., Kramer, C., Hitschfeld, M., Garcia-Burillo, S., & Mookerjea, B. 2007, A&A, 461, 143 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Sesto, L. A., Faifer, F. R., Smith Castelli, A. V., Forte, J. C., & Escudero, C. G. 2018, MNRAS, 479, 478 [NASA ADS] [Google Scholar]
Shibuya, T., Ouchi, M., Kubo, M., & Harikane, Y. 2016, ApJ, 821, 72 [NASA ADS] [CrossRef] [Google Scholar]
Silva-Villa, E., & Larsen, S. S. 2011, A&A, 529, A25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Somerville, R. S., & Davé, R. 2015, ARAA, 53, 51 [Google Scholar]
Somerville, R. S., Hopkins, P. F., Cox, T. J., Robertson, B. E., & Hernquist, L. 2008, MNRAS, 391, 481 [NASA ADS] [CrossRef] [Google Scholar]
Spengler, C., Côté, P., Roediger, J., et al. 2017, ApJ, 849, 55 [NASA ADS] [CrossRef] [Google Scholar]
Spina, L., Magrini, L., & Cunha, K. 2022, Universe, 8, 87 [NASA ADS] [CrossRef] [Google Scholar]
Spinoso, D., Bonoli, S., Valiante, R., Schneider, R., & Izquierdo-Villalba, D. 2023, MNRAS, 518, 4672 [Google Scholar]
Spitler, L. R., & Forbes, D. A. 2009, MNRAS, 392, L1 [Google Scholar]
Spitzer, L. 1940, MNRAS, 100, 396 [Google Scholar]
Spitzer, L. 1987, Dynamical Evolution of Globular Clusters (Princeton: Princeton University Press) [Google Scholar]
Springel, V., White, S. D. M., Tormen, G., & Kauffmann, G. 2001, MNRAS, 328, 726 [Google Scholar]
Springel, V., White, S. D. M., Jenkins, A., et al. 2005, Nature, 435, 629 [Google Scholar]
Strader, J., Brodie, J. P., Cenarro, A. J., Beasley, M. A., & Forbes, D. A. 2005, AJ, 130, 1315 [NASA ADS] [CrossRef] [Google Scholar]
Straižys, V., Maskoliūnas, M., Boyle, R. P., et al. 2014, MNRAS, 437, 1628 [Google Scholar]
Sun, J., Leroy, A. K., Rosolowsky, E., et al. 2022, AJ, 164, 43 [NASA ADS] [CrossRef] [Google Scholar]
Tanaka, T., & Haiman, Z. 2009, ApJ, 696, 1798 [NASA ADS] [CrossRef] [Google Scholar]
Thilker, D. A., Lee, J. C., Whitmore, B. C., et al. 2025, ApJSS, 280, 1 [Google Scholar]
Tonini, C. 2013, ApJ, 762, 39 [NASA ADS] [CrossRef] [Google Scholar]
Toomre, A. 1964, ApJ, 139, 1217 [Google Scholar]
Torrey, P., Vogelsberger, M., Marinacci, F., et al. 2019, MNRAS, 484, 5587 [NASA ADS] [Google Scholar]
Tremonti, C. A., Heckman, T. M., Kauffmann, G., et al. 2004, ApJ, 613, 898 [Google Scholar]
Trujillo-Gomez, S., Kruijssen, J. M. D., Reina-Campos, M., et al. 2021, MNRAS, 503, 31 [NASA ADS] [CrossRef] [Google Scholar]
Untzaga, J., Bonoli, S., Izquierdo-Villalba, D., Mezcua, M., & Spinoso, D. 2024, MNRAS, 535, 3293 [Google Scholar]
Usher, C., Forbes, D. A., Brodie, J. P., et al. 2012, MNRAS, 426, 1475 [Google Scholar]
Usher, C., Forbes, D. A., Spitler, L. R., et al. 2013, MNRAS, 436, 1172 [Google Scholar]
Usher, C., Caldwell, N., & Cabrera-Ziri, I. 2024, MNRAS, 528, 6010 [NASA ADS] [CrossRef] [Google Scholar]
Valenzuela, L. M., Moster, B. P., Remus, R.-S., O’Leary, J. A., & Burkert, A. 2021, MNRAS, 505, 5815 [NASA ADS] [CrossRef] [Google Scholar]
Valenzuela, L. M., Forbes, D. A., & Remus, R.-S. 2025, MNRAS, 537, 306 [CrossRef] [Google Scholar]
van der Kruit, P. C., & Freeman, K. C. 2011, ARAA, 49, 301 [Google Scholar]
Vani, A., Ayromlou, M., Kauffmann, G., & Springel, V. 2025, MNRAS, 536, 777 [NASA ADS] [Google Scholar]
Vanzella, E., Calura, F., Meneghetti, M., et al. 2017, MNRAS, 467, 4303 [Google Scholar]
Vanzella, E., Calura, F., Meneghetti, M., et al. 2019, MNRAS, 483, 3618 [Google Scholar]
Vanzella, E., Claeyssens, A., Welch, B., et al. 2023, ApJ, 945, 53 [NASA ADS] [CrossRef] [Google Scholar]
Vazquez-Semadeni, E. 1994, ApJ, 423, 681 [Google Scholar]
Vázquez-Semadeni, J. S. E., Chappell, D., & Passot, T. 1998, ApJ, 504, 835 [CrossRef] [Google Scholar]
Veljanoski, J., Ferguson, A. M. N., Mackey, A. D., et al. 2015, MNRAS, 452, 320 [NASA ADS] [CrossRef] [Google Scholar]
Vijayan, A. P., Scott, S. J., Thomas, P. A., et al. 2019, MNRAS, 489, 4072 [Google Scholar]
Villaume, A., Romanowsky, A. J., Brodie, J., & Strader, J. 2019, ApJ, 879, 45 [NASA ADS] [CrossRef] [Google Scholar]
Wang, K., & Peng, Y. 2024, A&A, 12 [Google Scholar]
Wang, E., Wang, H., Mo, H., et al. 2018, ApJ, 864, 51 [NASA ADS] [CrossRef] [Google Scholar]
Watson, G. N. 1944, A Treatise on the Theory of Bessel Functions, The Mathematical Gazette (Cambridge: Cambridge University Press) [Google Scholar]
West, M. J., Côté, P., Jones, C., Forman, W., & Marzke, R. O. 1995, ApJ, 453, L77 [NASA ADS] [CrossRef] [Google Scholar]
West, M. J., Côté, P., Marzke, R. O., & Jordán, A. 2004, Nature, 427, 31 [Google Scholar]
Westfall, K. B., Andersen, D. R., Bershady, M. A., et al. 2014, ApJ, 785, 43 [Google Scholar]
White, S. D. M., & Frenk, C. S. 1991, ApJ, 379, 52 [Google Scholar]
White, S. D. M., & Rees, M. J. 1978, MNRAS, 183, 341 [Google Scholar]
Whitmore, B. C., Chandar, R., Bowers, A. S., et al. 2014, AJ, 147, 78 [NASA ADS] [CrossRef] [Google Scholar]
Yates, R. M., Henriques, B. M. B., Thomas, P. A., et al. 2013, MNRAS, 435, 3500 [Google Scholar]
Yates, R. M., Thomas, P. A., & Henriques, B. M. B. 2017, MNRAS, 464, 3169 [NASA ADS] [CrossRef] [Google Scholar]
Yates, R. M., Henriques, B. M. B., Fu, J., et al. 2021, MNRAS, 503, 4474 [NASA ADS] [CrossRef] [Google Scholar]
Yates, R. M., Hendriks, D., Vijayan, A. P., et al. 2024, MNRAS, 527, 6292 [Google Scholar]
Zaritsky, D. 2022, MNRAS, 513, 2609 [NASA ADS] [CrossRef] [Google Scholar]
Zentner, A. R., & Bullock, J. S. 2003, ApJ, 598, 49 [NASA ADS] [CrossRef] [Google Scholar]
Zhang, Q., & Fall, S. M. 1999, ApJ, 527, L81 [CrossRef] [Google Scholar]
Zhong, W., Fu, J., Sharma, P., Shen, S., & Yates, R. M. 2023a, MNRAS, 519, 4344 [Google Scholar]
Zhong, W., Fu, J., Shen, S., & Yuan, F. 2023b, RAA, 23, 075004 [Google Scholar]
Zhou, L., Federrath, C., Yuan, T., et al. 2017, MNRAS, 470, 4573 [Google Scholar]

Yates et al. (2017, 2024) discuss and utilise more physically motivated profiles within L-Galaxies.

Using the same cluster formation efficiency calculations from Kruijssen (2012), Pfeffer et al. (2025) recently showed that the E-MOSAICS simulation can well reproduce observed properties of star clusters at high-z.

Here we completely neglect the increase in velocity dispersion as a function of stellar age caused by interactions with giant molecular clouds (GMCs) or spiral waves, potentially resulting in a more stable disk.

⁴

Note that the final distance, once τ_fric = ∆t typically is greater than zero (but often smaller than 1 kpc) due to the invalidity of the assumption of circular motion of the satellite around the central halo.

⁵

L-Galaxies does not model any spiral-wave patterns in galaxies, which is why we prefer the term “disk-dominated” over “spiral”.

⁶

Note that the secondary dependence on galaxy mass is also present in the E-MOSAICS model (Figure 10 in Pfeffer et al. 2019b); however, they only include 153 galaxies and have an additional y-axis offset of 0.5–1 dex.

⁷

The same trend is apparent for star clusters in bulge-dominated galaxies given their rich merger history compared to their disk-dominated counterparts.

⁸

There still remains some difference in the metallicities of newly formed star clusters between galaxies of different morphology, as we show in Section 3.5.

⁹

The dataset of Brown & Gnedin (2021) features only a few star clusters older than 6 Gyr. We omit a comparison for old star clusters.

¹⁰

Young star clusters in a galaxy’s halo do not form there. Instead, they formed in a satellite galaxy that was accreted shortly after their formation.

¹¹

However, the detection of a low-metallicity stellar stream of a former massive star cluster in the Milky Way halo (Martin et al. 2022) and the detection of a massive star cluster with [Fe/H] ≈ −2.9 in M 31 (Larsen et al. 2020) challenge this notion.

¹²

We confirmed that the random sampling of initial values during bootstrapping affects the bi-modal fractions at most by 0.05.

¹³

Note that we do not utilise a cut on the kurtosis of the Gaussians as we find in some cases that, if there is one dominant and another sub-dominant star cluster population in a galaxy, the kurtosis is greater than zero, indicating a ‘peaked’ distribution, thus rejecting the hypothesis of a bimodal distribution.

¹⁴

The parameter ϵ_c refers to the efficiency of star formation within a GMC and not the average efficiency of star formation within the interstellar medium that was introduced in Section 2.3.3.

Appendix A Model variations

We present here additional variations of model parameters that influence the assembly history of z = 0 star cluster populations.

Appendix A.1 Velocity dispersion of the cold gas

One of the most crucial parameters for modelling the Toomre parameter is the velocity dispersion of the cold gas. As discussed in e.g. Lehnert et al. (2009) and Zhou et al. (2017) there exists substantial scatter in the relationship between the velocity dispersion of the cold gas and the SFR surface density (c.f. Equation 14). In our fiducial model we adopted ***(74)*** $α_{cold}^{fid} = 5 {km s}^{- 1}$ $Mathematical equation: $\alpha _{{\rm{cold}}}^{{\rm{fid}}} = 5{\rm{km }}{{\rm{s}}^{ - 1}}$$ , ***(75)*** $β_{cold}^{fid} = 20 {km s}^{- 1}$ $Mathematical equation: $\beta _{{\rm{cold}}}^{{\rm{fid}}} = 20{\rm{km }}{{\rm{s}}^{ - 1}}$$ , and ***(76)*** $γ_{cold}^{fid} = 1 / 3$ $Mathematical equation: $\eqalign{& \beta _{{\rm{cold}}}^{{\rm{var}}} = 100{\rm{km }}{{\rm{s}}^{ - 1}} \cr} $$ for the offset, slope, and exponent, respectively.

Here we explore a wider parameter range by testing a first version where we increase the slope value to ***(77)*** $β_{cold}^{var} = 100 {km s}^{- 1}$ $Mathematical equation: $\beta _{{\rm{cold}}}^{{\rm{var}}} = 100{\rm{km }}{{\rm{s}}^{ - 1}}$$ and a second version where we set the exponent to ***(78)*** $γ_{cold}^{var} = 1 / 2$ $Mathematical equation: $\gamma _{{\rm{cold}}}^{{\rm{var}}} = 1/2$$ . Additionally, we test a separate prescription based on the Jeans mass (see e.g. Elmegreen et al. 2007) where $σ_{g, j} \sim M_{J}^{1 / 4} G^{1 / 2} Σ_{g, j}^{1 / 4} = 4.4 {pc Myr}^{- 1} Σ_{SFR, j}^{0.18},$ $Mathematical equation: ${\sigma _{g,j}} \sim M_{\rm{J}}^{1/4}{G^{1/2}}\Sigma _{g,j}^{1/4} = 4.4{\rm{pcMy}}{{\rm{r}}^{ - 1}}\Sigma _{{\rm{SFR,}}j}^{0.18},$$ (A.1)

with M_J as the Jeans mass, which we assume to be 10⁹ M_⊙ for the equality, and converted from the cold gas surface density to the SFR density using Equation 7 from Kennicutt (1998) with Σ_SFR in units of 10⁶ M_⊙pc⁻²Myr⁻¹. In addition, we add a velocity floor of 5 kms⁻¹, the same value that we used in Equation 14.

To probe the effect of this relationship on the properties of newly formed star clusters, we fit linear relationships to the resulting M_V-SFR parameter space, as presented in Section 3.1, and present the slope values in Table A.1.

We find almost no difference in the results between the fiducial model and the model where γ_cold = 1/2, indicating that the velocity dispersion of the cold gas has a subdominant effect on the effective Toomre stability parameter in comparison to the velocity dispersion of stars in the disk. For the case where we increase the importance of star formation by increasing the slope parameter β_cold we find a slightly steeper slope than in the fiducial model. This result highlights that the assumption on the relative importance of the star formation surface densities become increasingly important with lower-mass galaxies. Here the Toomre stability parameter of the cold gas disk gains more relative weight compared to the one of the stellar disk. Since Q_D,g ∝ c_s,cold, the Toomre parameter of the gas increases/decreases, which then also results in an increase/decrease of Q_eff, albeit not as strong as Q_D,g due to taking a weighted average with Q_D,s. This increase/decrease has a direct consequence on the star cluster masses because the upper truncation mass in Equations 28 and 29 scales as ***(80)*** $m_{cl,max} \propto Q_{eff}^{4}$ $Mathematical equation: ${m_{{\rm{cl,max}}}} \propto Q_{{\rm{eff}}}^4$$ , resulting in an increased/decreased probability to randomly sample massive clusters. The increase in stability is also expected given the mass dependence of the M_V-SFR relation that we present in Figure 4, which comes from the SFR versus galaxy mass relationship. For the same reason, we find a stark decline in the steepness of the relationship when adopting Equation A.1.

Our simulations indicate that the slope of the relationship is steeper for more massive galaxies irrespective of the assumed cold gas velocity dispersion versus surface SFR relationship.

Appendix A.2 Initial half-mass radius

The distribution of initial half-mass radii remains an unsolved problem. Measurements of the half-light size of young star clusters give typical values in the range of a few parsecs to about 0.5 pc (e.g. Bastian et al. 2012; Ryon et al. 2015; Brown & Gnedin 2021) but the scatter at similar cluster ages is large, going up to ≈ 20 pc for W 3 in NGC 7252 (Maraston et al. 2004; Fellhauer & Kroupa 2005). This could indicate that the initial half-mass radius of clusters is similar and diverges due to a different initial evolutionary phase given the local environment, or that other parameter influence the initial half-mass radius, such as cluster mass or the local physical conditions. Reina-Campos et al. (2023a) tested different prescriptions of initial half-mass radii, considering constant values, constant densities, a linear relationships from the data provided by Brown & Gnedin (2021), and a theoretical model from Choksi & Kruijssen (2021) that reads $r_{c, h} = {(\frac{3}{10 π^{2}} \frac{α_{vir}}{ϕ_{P} ϕ_{\bar{P}, j}} \frac{m_{c}^{2}}{Σ_{g}^{2}})}^{1 / 4} \frac{ε_{c}^{1 / 2}}{2 ε_{c} - 1} \frac{f_{acc}}{Q_{eff}^{2}},$ $Mathematical equation: ${r_{{\rm{c,}}h}} = {\left( {{3 \over {10{\pi ^2}}}{{{\alpha _{{\rm{vir}}}}} \over {{\phi _P}{\phi _{\bar P,j}}}}{{m_{\rm{c}}^2} \over {\Sigma _{\rm{g}}^2}}} \right)^{1/4}}{{\varepsilon _{\rm{c}}^{1/2}} \over {2{\varepsilon _c} - 1}}{{{f_{{\rm{acc}}}}} \over {Q_{{\rm{eff}}}^2}},$$ (A.2)

with f_acc = 0.6, ϵ_c = 1.0, and environmentally dependent parameters evaluated for each ring.¹⁴ Reina-Campos et al. (2023a) conclude that none of the prescriptions reproduce the size-mass relationship after evolving the cluster population for a few gigayears within the EMP-Pathfinder simulation suite (Reina-Campos et al. 2022a).

Table A.1

Slope values of the M_V-SFR relationship when adjusting parameters for the c_s,cold-Σ_SFR relationship.

Similar to Reina-Campos et al. (2023a) we implement different prescriptions for the initial half-mass radius, that we outline in Table A.2, to test their influence on the distributions of stellar mass and half-mass radii, as discussed in Section 3.2 and 3.3.

We find that the precise details of the initial half-mass radius affects the star cluster properties (see Figure A.1 and A.2 below). When we choose a more compact initial half-mass radius, the star clusters expand more quickly at young ages, roughly matching the distribution of young star clusters in the observational dataset of Brown & Gnedin (2021) in dwarf galaxies. However the agreement in the half-mass radii becomes worse in more massive galaxies and the slope of the M_v-SFR relationship steepens as less massive star clusters loose more mass, resulting in an increase of M_V at low SFRs.

When using Eq. A.2, we find a similar distribution compared to the fiducial model as many star clusters have initial half-mass radii only slightly smaller than 0.1 pc given the weak dependence on star cluster mass and the inverse squared dependence on the Toomre stability parameter. As a consequence, similar to the fiducial model, this choice of initial half-mass radii fails to reproduce the observational constraints by Brown & Gnedin (2021).

Table A.2

Variations in the initial half-mass radii of star clusters.

Looking at the masses of the star clusters we find stronger deviations for the “compact” and “theoretical” scenario. In the first one, star clusters start out compact and expand quickly in their dense environment, which results in an increase in mass loss. In turn, this effect shifts the distribution of star cluster masses towards lower values. For the second scenario, the initial half-mass radius is proportional to the square-root of the star cluster mass. Therefore, low-mass star clusters expand rapidly and lose mass quickly, directly affecting the mass distribution. Note that, as the environment is denser in more massive disks, the mass distribution is most strongly affected in massive disk galaxies.

Appendix A.3 Cluster initial mass function

To test the effect of the upper truncation mass on the resulting star cluster populations we re-run the model considering a pure power-law CIMF, i.e. ***(94)*** $ξ (m_{c}) \propto m_{c}^{α}$ $Mathematical equation: $\xi \left( {{m_c}} \right) \propto m_{\rm{c}}^\alpha $$ .

While the overall shape of the half-mass radius distribution remains unchanged, we find that star cluster masses old the old population are elevated by 0.2 to 0.3 dex. This rather small difference of a factor of two in star cluster mass comes from the relatively high upper truncation mass of the CIMF at high-z. A small difference of the same magnitude is still present in the population of young star clusters. However, the distribution has a longer tail towards higher stellar masses compared to the fiducial model.

This result translated directly into the M_V-SFR relationship as well where the slope value increases (from α ≈ −2.47 for the fiducial model to α ≈ −2.33 for the model variation) and absolute magnitudes are slightly increased. As a consequence, the model variation now fails to account for some of the lowest-mass star clusters in this parameter space around SFRs of 1 M_⊙ yr⁻¹. Other parameter distributions, such as the mean star cluster metallicity versus galaxy mass remain unchanged.

Fig. A.1

Comparison between the distribution of half-mass radii of young star clusters in the disks of disk-dominated galaxies (see Figure 5 for details) when considering different assumptions on the initial values, as detailed in Table A.2. Black solid lines give the observational dataset of Brown & Gnedin (2021).

Fig. A.2

Same as Figure A.1 but for star cluster masses.

All Tables

Table A.1

Slope values of the M_V-SFR relationship when adjusting parameters for the c_s,cold-Σ_SFR relationship.

In the text

Table A.2

Variations in the initial half-mass radii of star clusters.

In the text

All Figures

	Fig. 1 Relevant prescriptions for the assembly and evolution of star clusters implemented into a modified version (Yates et al. 2021) of L-Galaxies 2020 (Henriques et al. 2020).
In the text

	Fig. 6 Same as Figure 5 but for star cluster half-mass radii. The double-peaked structure of the half-mass radii in a galaxy’s halo is a result of the prescription for star cluster re-distribution during galaxy mergers (see Sections 2.5.4 and 3.4 below).
In the text

	Fig. 7 Same as Figure 5 but for star cluster half-mass radii. The double-peaked structure of the half-mass radii in a galaxy’s halo is a result of the prescription for star cluster re-distribution during galaxy mergers.
In the text

	Fig. A.1 Comparison between the distribution of half-mass radii of young star clusters in the disks of disk-dominated galaxies (see Figure 5 for details) when considering different assumptions on the initial values, as detailed in Table A.2. Black solid lines give the observational dataset of Brown & Gnedin (2021).
In the text

	Fig. A.2 Same as Figure A.1 but for star cluster masses.
In the text

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

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

Initial download of the metrics may take a while.

[R1] Adamo, A., Östlin, G., & Zackrisson, E. 2011, MNRAS, 417, 1904 [NASA ADS] [CrossRef] [Google Scholar]

[R2] Adamo, A., Kruijssen, J. M. D., Bastian, N., Silva-Villa, E., & Ryon, J. 2015, MNRAS, 452, 246 [NASA ADS] [CrossRef] [Google Scholar]

[R3] Adamo, A., Ryon, J. E., Messa, M., et al. 2017, ApJ, 841, 131 [Google Scholar]

[R4] Adamo, A., Bradley, L. D., Vanzella, E., et al. 2024, Nature, 632, 513 [NASA ADS] [CrossRef] [Google Scholar]

[R5] Alexander, P. E. R., & Gieles, M. 2012, MNRAS, 422, 3415 [NASA ADS] [CrossRef] [Google Scholar]

[R6] Alexander, P. E. R., Gieles, M., Lamers, H. J. G. L. M., & Baumgardt, H. 2014, MNRAS, 442, 1265 [Google Scholar]

[R7] Andersson, E. P., Mac Low, M.-M., Agertz, O., Renaud, F., & Li, H. 2024, A&A, 681, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[R8] Angulo, R. E., & Hilbert, S. 2015, MNRAS, 448, 364 [CrossRef] [Google Scholar]

[R9] Angulo, R. E., & White, S. D. M. 2010, MNRAS, 405, 143 [NASA ADS] [Google Scholar]

[R10] Annibali, F., Tosi, M., Monelli, M., et al. 2009, AJ, 138, 169 [NASA ADS] [CrossRef] [Google Scholar]

[R11] Annibali, F., Tosi, M., Aloisi, A., & van der Marel, R. P. 2011, AJ, 142, 129 [NASA ADS] [CrossRef] [Google Scholar]

[R12] Ayromlou, M., Kauffmann, G., Yates, R. M., Nelson, D., & White, S. D. M. 2021, MNRAS, 505, 492 [NASA ADS] [CrossRef] [Google Scholar]

[R13] Bajkova, A. T., Smirnov, A. A., & Bobylev, V. V. 2023, AB, 78, 499 [Google Scholar]

[R14] Banerjee, S. 2017, MNRAS, 467, 524 [NASA ADS] [Google Scholar]

[R15] Banerjee, S., & Kroupa, P. 2017, A&A, 597, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[R16] Barmby, P., & Huchra, J. P. 2001, AJ, 122, 2458 [CrossRef] [Google Scholar]

[R17] Barmby, P., Huchra, J. P., Brodie, J. P., et al. 2000, AJ, 119, 727 [NASA ADS] [CrossRef] [Google Scholar]

[R18] Barrera, M., Springel, V., White, S. D. M., et al. 2023, MNRAS, 525, 6312 [NASA ADS] [CrossRef] [Google Scholar]

[R19] Bassino, L. P., & Caso, J. P. 2017, MNRAS, 466, 4259 [NASA ADS] [Google Scholar]

[R20] Bastian, N. 2008, MNRAS, 390, 759 [Google Scholar]

[R21] Bastian, N., Adamo, A., Gieles, M., et al. 2012, MNRAS, 419, 2606 [Google Scholar]

[R22] Bastian, N., Pfeffer, J., Kruijssen, J. M. D., et al. 2020, MNRAS, 498, 1050 [Google Scholar]

[R23] Baugh, C. M. 2006, Rep. Prog. Phys., 69, 3101 [NASA ADS] [CrossRef] [Google Scholar]

[R24] Baumgardt, H., & Hilker, M. 2018, MNRAS, 478, 1520 [Google Scholar]

[R25] Beasley, M. A., Bridges, T., Peng, E., et al. 2008, MNRAS, 386, 1443 [NASA ADS] [CrossRef] [Google Scholar]

[R26] Beasley, M. A., Leaman, R., Gallart, C., et al. 2019, MNRAS, 487, 1986 [Google Scholar]

[R27] Bertin, G., & Romeo, A. B. 1988, A&A, 195, 105 [NASA ADS] [Google Scholar]

[R28] Bertoldi, F., & McKee, C. F. 1992, ApJ, 395, 140 [NASA ADS] [CrossRef] [Google Scholar]

[R29] Bik, A., Lamers, H. J. G. L. M., Bastian, N., Panagia, N., & Romaniello, M. 2003, A&A, 397, 473 [CrossRef] [EDP Sciences] [Google Scholar]

[R30] Binney, J., & Tremaine, S. 2008, Galactic Dynamics, 2nd edn. (Princeton: Princeton University Press) [Google Scholar]

[R31] Blakeslee, J. P. 1999, AJ, 118, 1506 [NASA ADS] [CrossRef] [Google Scholar]

[R32] Blom, C., Spitler, L. R., & Forbes, D. A. 2012, MNRAS, 420, 37 [Google Scholar]

[R33] Bluck, A. F. L., Mendel, T. J., Ellison, S. L., et al. 2016, MNRAS, 462, 2559 [Google Scholar]

[R34] Bonoli, S., Marulli, F., Springel, V., et al. 2009, MNRAS, 396, 423 [NASA ADS] [CrossRef] [Google Scholar]

[R35] Bonoli, S., Mayer, L., & Callegari, S. 2014, MNRAS, 437, 1576 [NASA ADS] [CrossRef] [Google Scholar]

[R36] Bottema, R. 1993, A&A, 275, 16 [NASA ADS] [Google Scholar]

[R37] Boylan-Kolchin, M., Springel, V., White, S. D. M., Jenkins, A., & Lemson, G. 2009, MNRAS, 398, 1150 [Google Scholar]

[R38] Bradley, L. D., Adamo, A., Vanzella, E., et al. 2024, ApJ, 22 [Google Scholar]

[R39] Brodie, J. P., & Strader, J. 2006, ARAA, 44, 193 [Google Scholar]

[R40] Brodie, J. P., Usher, C., Conroy, C., et al. 2012, ApJ, 759, L33 [NASA ADS] [CrossRef] [Google Scholar]

[R41] Brogaard, K., VandenBerg, D. A., Bruntt, H., et al. 2012, A&A, 543, A106 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[R42] Brown, G., & Gnedin, O. Y. 2021, MNRAS, 508, 5935 [NASA ADS] [CrossRef] [Google Scholar]

[R43] Brown, G., & Gnedin, O. Y. 2022, MNRAS, 514, 280 [Google Scholar]

[R44] Burkert, A., & Forbes, D. A. 2020, AJ, 159, 56 [Google Scholar]

[R45] Burkhart, B. 2018, ApJ, 863, 118 [CrossRef] [Google Scholar]

[R46] Caldwell, N., & Romanowsky, A. J. 2016, ApJ, 824, 42 [NASA ADS] [CrossRef] [Google Scholar]

[R47] Callingham, T. M., Cautun, M., Deason, A. J., et al. 2022, MNRAS, 513, 4107 [NASA ADS] [CrossRef] [Google Scholar]

[R48] Calura, F., Pascale, R., Agertz, O., et al. 2025, A&A, 698, A207 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[R49] Calzetti, D., Lee, J. C., Sabbi, E., et al. 2015, AJ, 149, 51 [Google Scholar]

[R50] Cantiello, M., Capaccioli, M., Napolitano, N., et al. 2015, A&A, 576, A14 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[R51] Carretta, E., Gratton, R. G., Clementini, G., & Fusi Pecci, F. 2000, ApJ, 533, 215 [NASA ADS] [CrossRef] [Google Scholar]

[R52] Chabrier, G. 2003, PASP, 115, 763 [Google Scholar]

[R53] Chandrasekhar, S. 1943, ApJ, 97, 255 [Google Scholar]

[R54] Chattopadhyay, D., Hurley, J., Stevenson, S., & Raidani, A. 2022, MNRAS, 513, 4527 [NASA ADS] [CrossRef] [Google Scholar]

[R55] Chen, Y., Mo, H., & Wang, H. 2025, MNRAS, 540, 1235 [Google Scholar]

[R56] Chevance, M., Kruijssen, J. M. D., Hygate, A. P. S., et al. 2020, MNRAS, 493, 2872 [NASA ADS] [CrossRef] [Google Scholar]

[R57] Choksi, N., & Gnedin, O. Y. 2019, MNRAS, 486, 331 [NASA ADS] [CrossRef] [Google Scholar]

[R58] Choksi, N., & Kruijssen, J. M. D. 2021, MNRAS, 507, 5492 [Google Scholar]

[R59] Choksi, N., Behroozi, P., Volonteri, M., et al. 2017, MNRAS, 472, 1526 [Google Scholar]

[R60] Claeyssens, A., Adamo, A., Richard, J., et al. 2023, MNRAS, 520, 2180 [NASA ADS] [CrossRef] [Google Scholar]

[R61] Cohen, J. G., Blakeslee, J. P., & Ryzhov, A. 1998, ApJ, 496, 808 [NASA ADS] [CrossRef] [Google Scholar]

[R62] Colombo, D., Hughes, A., Schinnerer, E., et al. 2014, ApJ, 784, 3 [NASA ADS] [CrossRef] [Google Scholar]

[R63] Conroy, C., & Gunn, J. E. 2010a, Astrophysics Source Code Library [record ascl:1010.043] [Google Scholar]

[R64] Conroy, C., & Gunn, J. E. 2010b, ApJ, 712, 833 [Google Scholar]

[R65] Conroy, C., White, M., & Gunn, J. E. 2010, ApJ, 708, 58 [NASA ADS] [CrossRef] [Google Scholar]

[R66] Conroy, C., Gunn, J. E., & White, M. 2009, ApJ, 699, 486 [Google Scholar]

[R67] Conselice, C. J., Grogin, N. A., Jogee, S., et al. 2004, ApJ, 600, L139 [NASA ADS] [CrossRef] [Google Scholar]

[R68] Cook, D. O., Seth, A. C., Dale, D. A., et al. 2012, ApJ, 751, 100 [NASA ADS] [CrossRef] [Google Scholar]

[R69] Cook, D. O., Lee, J. C., Adamo, A., et al. 2023, MNRAS, 519, 3749 [CrossRef] [Google Scholar]

[R70] Correnti, M., Gennaro, M., Kalirai, J. S., Brown, T. M., & Calamida, A. 2016, ApJ, 823, 18 [Google Scholar]

[R71] Crain, R. A., Schaye, J., Bower, R. G., et al. 2015, MNRAS, 450, 1937 [NASA ADS] [CrossRef] [Google Scholar]

[R72] Croton, D. J. 2006, MNRAS, 369, 1808 [NASA ADS] [CrossRef] [Google Scholar]

[R73] Croton, D. J., Springel, V., White, S. D. M., et al. 2006, MNRAS, 365, 11 [Google Scholar]

[R74] da Silva, R. L., Fumagalli, M., & Krumholz, M. 2012, ApJ, 745, 145 [CrossRef] [Google Scholar]

[R75] da Silva, R. L., Fumagalli, M., & Krumholz, M. 2014, MNRAS, 444, 3275 [Google Scholar]

[R76] Davies, B., de La Fuente, D., Najarro, F., et al. 2012, MNRAS, 419, 1860 [Google Scholar]

[R77] Davis, M., Efstathiou, G., Frenk, C. S., & White, S. D. M. 1985, ApJ, 292, 371 [Google Scholar]

[R78] De Lucia, G., & Blaizot, J. 2007, MNRAS, 375, 2 [Google Scholar]

[R79] De Lucia, G., Tornatore, L., Frenk, C. S., et al. 2014, MNRAS, 445, 970 [Google Scholar]

[R80] De Lucia, G., Kruijssen, J. M. D., Trujillo-Gomez, S., Hirschmann, M., & Xie, L. 2024, MNRAS, 530, 2760 [NASA ADS] [CrossRef] [Google Scholar]

[R81] del P. Lagos, C., Bravo, M., Tobar, R., et al. 2024, MNRAS, 531, 3551 [Google Scholar]

[R82] Dillamore, A. M., Monty, S., Belokurov, V., & Evans, N. W. 2024, ApJ, 971, L4 [Google Scholar]

[R83] Dobbs, C. L., Burkert, A., & Pringle, J. E. 2011, MNRAS, 417, 1318 [NASA ADS] [CrossRef] [Google Scholar]

[R84] Dolag, K., Borgani, S., Murante, G., & Springel, V. 2009, MNRAS, 399, 497 [Google Scholar]

[R85] Dolfi, A., Pfeffer, J., Forbes, D. A., et al. 2022, MNRAS, 511, 3179 [Google Scholar]

[R86] Downing, J. M. B. 2012, MNRAS, 425, 2234 [Google Scholar]

[R87] Dutton, A. A., & Macció, A. V. 2014, MNRAS, 441, 3359 [NASA ADS] [CrossRef] [Google Scholar]

[R88] Ekström, S., Georgy, C., Eggenberger, P., et al. 2012, A&A, 537, A146 [Google Scholar]

[R89] El-Badry, K., Quataert, E., Weisz, D. R., Choksi, N., & Boylan-Kolchin, M. 2019, MNRAS, 482, 4528 [Google Scholar]

[R90] Elmegreen, B., & Lahén, N. 2024, OJA, 7, 68 [Google Scholar]

[R91] Elmegreen, D. M., Elmegreen, B. G., Ravindranath, S., & Coe, D. A. 2007, ApJ, 658, 763 [Google Scholar]

[R92] Emig, K. L., Bolatto, A. D., Leroy, A. K., et al. 2020, ApJ, 903, 50 [NASA ADS] [CrossRef] [Google Scholar]

[R93] Emsellem, E., Schinnerer, E., Santoro, F., et al. 2022, A&A, 659, A191 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[R94] Escala, A., Larson, R. B., Coppi, P. S., & Mardones, D. 2004, ApJ, 607, 765 [NASA ADS] [CrossRef] [Google Scholar]

[R95] Escudero, C. G., Faifer, F. R., Bassino, L. P., Calderón, J. P., & Caso, J. P. 2015, MNRAS, 449, 612 [Google Scholar]

[R96] Fahrion, K., Lyubenova, M., Hilker, M., et al. 2020, A&A, 637, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[R97] Fahrion, K., Lyubenova, M., van de Ven, G., et al. 2021, A&A, 650, A137 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[R98] Fellhauer, M., & Kroupa, P. 2005, MNRAS, 359, 223 [Google Scholar]

[R99] Feng, C.-C., Lin, L.-H., Wang, H.-H., & Taam, R. E. 2014, ApJ, 785, 103 [Google Scholar]

[R100] Filion, C., McClure, R. L., Weinberg, M. D., D’Onghia, E., & Daniel, K. J. 2023, MNRAS, 524, 276 [Google Scholar]

[R101] Flynn, C., Holmberg, J., Portinari, L., Fuchs, B., & Jahreiß, H. 2006, MNRAS, 372, 1149 [NASA ADS] [CrossRef] [Google Scholar]

[R102] Forbes, D. A., & Bridges, T. 2010, MNRAS, 404, 1203 [NASA ADS] [Google Scholar]

[R103] Forbes, D. A., & Romanowsky, A. J. 2023, MNRAS, 520, L58 [Google Scholar]

[R104] Forbes, D. A., Bastian, N., Gieles, M., et al. 2018, Proc. R. Soc. London Ser. A, 474, 20170616 [Google Scholar]

[R105] Förster Schreiber, N. M., Shapley, A. E., Erb, D. K., et al. 2011, ApJ, 731, 65 [CrossRef] [Google Scholar]

[R106] Freeman, K. C. 1970, ApJ, 160, 811 [Google Scholar]

[R107] Fu, J., Guo, Q., Kauffmann, G., & Krumholz, M. R. 2010, MNRAS, 409, 515 [Google Scholar]

[R108] Fu, J., Kauffmann, G., Li, C., & Qi, G. 2012, MNRAS, 424, 2701 [Google Scholar]

[R109] Fu, J., Kauffmann, G., Huang, M.-L., et al. 2013, MNRAS, 434, 1531 [NASA ADS] [CrossRef] [Google Scholar]

[R110] Fukui, Y., Ohama, A., Hanaoka, N., et al. 2014, ApJ, 780, 36 [Google Scholar]

[R111] Gabrielpillai, A., Somerville, R. S., Genel, S., et al. 2022, MNRAS, 517, 6091 [Google Scholar]

[R112] Garcia, F. A. B., Ricotti, M., Sugimura, K., & Park, J. 2023, MNRAS, 522, 2495 [NASA ADS] [CrossRef] [Google Scholar]

[R113] Garro, E. R., Minniti, D., & Fernández-Trincado, J. G. 2024, A&A, 687, A214 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[R114] Genzel, R., Newman, S., Jones, T., et al. 2011, ApJ, 733, 101 [Google Scholar]

[R115] Georgy, C., Ekström, S., Eggenberger, P., et al. 2013, A&A, 558, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[R116] Gieles, M., & Renaud, F. 2016, MNRAS, 463, L103 [Google Scholar]

[R117] Gieles, M., Heggie, D. C., & Zhao, H. 2011, MNRAS, 413, 2509 [Google Scholar]

[R118] Gieles, M., Erkal, D., Antonini, F., Balbinot, E., & Peñarrubia, J. 2021, NatAs, 5, 957 [Google Scholar]

[R119] Giersz, M., & Heggie, D. C. 1994, MNRAS, 268, 257 [NASA ADS] [CrossRef] [Google Scholar]

[R120] Giersz, M., Askar, A., Wang, L., et al. 2019, MNRAS, 487, 2412 [NASA ADS] [CrossRef] [Google Scholar]

[R121] Goddard, Q. E., Bastian, N., & Kennicutt, R. C. 2010, MNRAS, 405, 857 [NASA ADS] [Google Scholar]

[R122] Green, A. W., Glazebrook, K., McGregor, P. J., et al. 2014, MNRAS, 437, 1070 [Google Scholar]

[R123] Guo, S., & Qi, Z. 2023, Universe, 9, 252 [Google Scholar]

[R124] Guo, Q., White, S., Boylan-Kolchin, M., et al. 2011, MNRAS, 413, 101 [Google Scholar]

[R125] Gustafsson, B., Church, R. P., Davies, M. B., & Rickman, H. 2016, A&A, 593, A85 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[R126] Hackshaw, Z., Hawkins, K., Filion, C., et al. 2024, ApJ, 255, 19 [Google Scholar]

[R127] Hammer, F., Li, H., Mamon, G. A., et al. 2023, MNRAS, 519, 5059 [CrossRef] [Google Scholar]

[R128] Hammer, F., Wang, J., Mamon, G. A., et al. 2024, MNRAS, 527, 2718 [Google Scholar]

[R129] Harris, W. E., Ciccone, S. M., Eadie, G. M., et al. 2017, ApJ, 835, 101 [NASA ADS] [CrossRef] [Google Scholar]

[R130] Hawkins, K. 2023, MNRAS, 525, 3318 [NASA ADS] [CrossRef] [Google Scholar]

[R131] Heiter, U., Soubiran, C., Netopil, M., & Paunzen, E. 2014, A&A, 561, A93 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[R132] Hénon, M. 1961, Annales d’Astrophysique, 24, 369 [NASA ADS] [Google Scholar]

[R133] Henriques, B. M. B., White, S. D. M., Thomas, P. A., et al. 2015, MNRAS, 451, 2663 [Google Scholar]

[R134] Henriques, B. M. B., Yates, R. M., Fu, J., et al. 2020, MNRAS, 491, 5795 [NASA ADS] [CrossRef] [Google Scholar]

[R135] Heyer, M. H., Corbelli, E., Schneider, S. E., & Young, J. S. 2004, ApJ, 602, 723 [NASA ADS] [CrossRef] [Google Scholar]

[R136] Heyer, M., Krawczyk, C., Duval, J., & Jackson, J. M. 2009, ApJ, 699, 1092 [Google Scholar]

[R137] Hirschmann, M., De Lucia, G., & Fontanot, F. 2016, MNRAS, 461, 1760 [Google Scholar]

[R138] Hislop, J. M., Naab, T., Steinwandel, U. P., et al. 2022, MNRAS, 509, 5938 [Google Scholar]

[R139] Hixenbaugh, K., Chandar, P., & Mok, A. 2022, AJ, 163, 271 [Google Scholar]

[R140] Ho, I.-T., Seibert, M., Meidt, S. E., et al. 2017, ApJ, 846, 39 [NASA ADS] [CrossRef] [Google Scholar]

[R141] Ho, I.-T., Kreckel, K., Meidt, S. E., et al. 2019, ApJ, 885, L31 [NASA ADS] [CrossRef] [Google Scholar]

[R142] Hollyhead, K., Adamo, A., Bastian, N., Gieles, M., & Ryon, J. E. 2016, MNRAS, 460, 2087 [Google Scholar]

[R143] Hopkins, P. F., Quataert, E., & Murray, N. 2012, MNRAS, 421, 3488 [NASA ADS] [CrossRef] [Google Scholar]

[R144] Horta, D., Hughes, M. E., Pfeffer, J. L., et al. 2021, MNRAS, 500, 4768 [Google Scholar]

[R145] Hoyer, N., Neumayer, N., Georgiev, I. Y., Seth, A. C., & Greene, J. E. 2021, MNRAS, 507, 3246 [NASA ADS] [CrossRef] [Google Scholar]

[R146] Huchra, J. P., Brodie, J. P., & Kent, S. M. 1991, ApJ, 370, 495 [Google Scholar]

[R147] Hudson, M. J., Harris, G. L., & Harris, W. E. 2014, ApJ, 787, L5 [Google Scholar]

[R148] Huertas-Company, M., Iyer, K. G., Angeloudi, E., et al. 2024, A&A, 685, A48 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[R149] Hughes, M. E., Pfeffer, J., Martig, M., et al. 2019, MNRAS, 482, 2795 [Google Scholar]

[R150] Hughes, M. E., Pfeffer, J. L., Martig, M., et al. 2020, MNRAS, 491, 4012 [NASA ADS] [CrossRef] [Google Scholar]

[R151] Hughes, M. E., Pfeffer, J. L., Bastian, N., et al. 2022, MNRAS, 510, 6190 [Google Scholar]

[R152] Hunt, E. L., & Reffert, S. 2023, A&A, 673, A114 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[R153] Hunter, D. A., Elmegreen, B. G., Dupuy, T. J., & Mortonson, M. 2003, AJ, 126, 1836 [NASA ADS] [CrossRef] [Google Scholar]

[R154] Huxor, A. P., Mackey, A. D., Ferguson, A. M. N., et al. 2014, MNRAS, 442, 2165 [Google Scholar]

[R155] Hwang, H.-C., Barrera-Ballesteros, J. K., Heckman, T. M., et al. 2019, ApJ, 872, 144 [NASA ADS] [CrossRef] [Google Scholar]

[R156] Ines Ennis, A., Caso, J. P., & Bassino, L. P. 2024, OJA, 12 [Google Scholar]

[R157] Irodotou, D., Thomas, P. A., Henriques, B. M. B., Sargent, M. T., & Hislop, J. M. 2019, MNRAS, 489, 3609 [Google Scholar]

[R158] Izquierdo-Villalba, D., Bonoli, S., Spinoso, D., et al. 2019, MNRAS, 488, 609 [NASA ADS] [CrossRef] [Google Scholar]

[R159] Izquierdo-Villalba, D., Bonoli, S., Dotti, M., et al. 2020, MNRAS, 495, 4681 [NASA ADS] [CrossRef] [Google Scholar]

[R160] Izquierdo-Villalba, D., Sesana, A., Bonoli, S., & Colpi, M. 2022, MNRAS, 509, 3488 [Google Scholar]

[R161] Jaffe, W. 1983, MNRAS, 202, 995 [NASA ADS] [Google Scholar]

[R162] Johnson, B., Foreman-Mackey, D., Sick, J., et al. 2023, https://doi.org/10.5281/zenodo.10026684 [Google Scholar]

[R163] Johnson, K. E., Leitherer, C., Vacca, W. D., & Conti, P. S. 2000, AJ, 120, 1273 [NASA ADS] [CrossRef] [Google Scholar]

[R164] Kacharov, N., Neumayer, N., Seth, A. C., et al. 2018, MNRAS, 480, 1973 [Google Scholar]

[R165] Karam, J., & Sills, A. 2022, MNRAS, 513, 6095 [Google Scholar]

[R166] Karam, J., & Sills, A. 2023, MNRAS, 521, 5557 [Google Scholar]

[R167] Katz, H., & Ricotti, M. 2013, MNRAS, 432, 3250 [Google Scholar]

[R168] Kauffmann, G., & Haehnelt, M. 2000, MNRAS, 311, 576 [Google Scholar]

[R169] Kaviraj, S., Ferreras, I., Yoon, S.-J., & Yi, S. K. 2005, A&A, 439, 913 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[R170] Keller, B. W., Kruijssen, J. M. D., Pfeffer, J., et al. 2020, MNRAS, 495, 4248 [NASA ADS] [CrossRef] [Google Scholar]

[R171] Kenney, J. D. P., & Lord, S. D. 1991, ApJ, 381, 118 [NASA ADS] [CrossRef] [Google Scholar]

[R172] Kennicutt, R. C. 1998, ARAA, 36, 189 [NASA ADS] [CrossRef] [Google Scholar]

[R173] Kewley, L. J., & Ellison, S. L. 2008, ApJ, 681, 1183 [Google Scholar]

[R174] Kikuchihara, S., Ouchi, M., Ono, Y., et al. 2020, ApJ, 893, 60 [NASA ADS] [CrossRef] [Google Scholar]

[R175] King, I. 1962, AJ, 67, 471 [Google Scholar]

[R176] Krauss, L. M., & Chaboyer, B. 2003, Science, 229, 65 [Google Scholar]

[R177] Kravtsov, A. V., & Gnedin, O. Y. 2005, ApJ, 623, 650 [Google Scholar]

[R178] Kreckel, K., Ho, I.-T., Blanc, G. A., et al. 2019, ApJ, 887, 80 [NASA ADS] [CrossRef] [Google Scholar]

[R179] Kritsuk, A. G., Norman, M. L., Padoan, P., & Wagner, R. 2007, ApJ, 665, 416 [NASA ADS] [CrossRef] [Google Scholar]

[R180] Kritsuk, A. G., Ustyugov, S. D., & Norman, M. L. 2017, New J. Phys., 19, 065003 [Google Scholar]

[R181] Kruijssen, J. M. D. 2012, MNRAS, 426, 3008 [Google Scholar]

[R182] Kruijssen, J. M. D. 2014, Class. Quantum Gravity, 31, 244006 [NASA ADS] [CrossRef] [Google Scholar]

[R183] Kruijssen, J. M. D. 2015, MNRAS, 454, 1658 [Google Scholar]

[R184] Kruijssen, J. M. D., & Cooper, A. P. 2012, MNRAS, 420, 340 [Google Scholar]

[R185] Kruijssen, J. M. D., Pfeffer, J. L., Crain, R. A., & Bastian, N. 2019a, MNRAS, 486, 3134 [NASA ADS] [CrossRef] [Google Scholar]

[R186] Kruijssen, J. M. D., Pfeffer, J. L., Reina-Campos, M., Crain, R. A., & Bastian, N. 2019b, MNRAS, 486, 3180 [Google Scholar]

[R187] Kruijssen, J. M. D., Pfeffer, J. L., Chevance, M., et al. 2020, MNRAS, 498, 2472 [NASA ADS] [CrossRef] [Google Scholar]

[R188] Krumholz, M. R., & Burkhart, B. 2016, MNRAS, 458, 1671 [NASA ADS] [CrossRef] [Google Scholar]

[R189] Krumholz, M. R., & McKee, C. F. 2005, ApJ, 630, 260 [Google Scholar]

[R190] Krumholz, M. R., McKee, C. F., & Tumlinson, J. 2009, ApJ, 693, 216 [NASA ADS] [CrossRef] [Google Scholar]

[R191] Krumholz, M. R., Fumagalli, M., da Silva, R. L., Rendahl, T., & Parra, J. 2015, MNRAS, 452, 1447 [NASA ADS] [CrossRef] [Google Scholar]

[R192] Krumholz, M. R., Burkhart, B., Forbes, J. C., & Crocker, R. M. 2018, MNRAS, 477, 2716 [CrossRef] [Google Scholar]

[R193] Kuijken, K., & Gilmore, G. 1989, MNRAS, 239, 571 [NASA ADS] [CrossRef] [Google Scholar]

[R194] Kundu, A., & Whitmore, B. C. 2001, AJ, 121, 2950 [NASA ADS] [CrossRef] [Google Scholar]

[R195] Lahén, N., Naab, T., Johansson, P. H., et al. 2020, ApJ, 891, 2 [CrossRef] [Google Scholar]

[R196] Lahén, N., Naab, T., Kauffmann, G., et al. 2023, MNRAS, 522, 3092 [CrossRef] [Google Scholar]

[R197] Lahén, N., Naab, T., & Szécsi, D. 2024, MNRAS, 530, 645 [CrossRef] [Google Scholar]

[R198] Lahén, N., Rantala, A., Naab, T., et al. 2025, MNRAS, 538, 2129 [Google Scholar]

[R199] Lake, W., Noaz, S., Chiou, Y. S., et al. 2021, ApJ, 922, 86 [Google Scholar]

[R200] Lake, W., Naoz, S., Marinacci, F., et al. 2023, ApJ, 956, L7 [Google Scholar]

[R201] Lamers, H. J. G. L. M., Gieles, M., & Portegies Zwart, S. F. 2005, A&A, 429, 173 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[R202] Larsen, S. S. 2002, AJ, 124, 1393 [NASA ADS] [CrossRef] [Google Scholar]

[R203] Larsen, S. S. 2010, Philos. Trans. R. Soc. A, 368, 867 [Google Scholar]

[R204] Larsen, S. S., Romanowsky, A. J., Brodie, J. P., & Wasserman, A. 2020, Science, 370, 970 [NASA ADS] [CrossRef] [Google Scholar]

[R205] Larson, R. B. 1981, MNRAS, 194, 809 [Google Scholar]

[R206] Le, M. N., & Cooper, A. P. 2025, ApJ, 978, 33 [Google Scholar]

[R207] Leaman, R., VandenBerg, D. A., & Mendel, J. T. 2013, MNRAS, 436, 122 [Google Scholar]

[R208] Lee, J. C., Sandstrom, K. M., Leroy, A. K., et al. 2023, ApJ, 944, L17 [NASA ADS] [CrossRef] [Google Scholar]

[R209] Lee, J. C., Whitmore, B. C., Thilker, D. A., et al. 2022, ApJSS, 258, 10 [Google Scholar]

[R210] Lee, M. G. 2003, JKAS, 36, 189 [Google Scholar]

[R211] Lehnert, M. D., Nesvadba, N. P. H., Le Tiran, L., et al. 2009, ApJ, 699, 1660 [NASA ADS] [CrossRef] [Google Scholar]

[R212] Leroy, A. K., Walter, F., Brinks, E., et al. 2008, AJ, 136, 2782 [Google Scholar]

[R213] Levy, R. C., Bolatto, A. D., Mayya, D., et al. 2024, ApJ, 973, 16 [Google Scholar]

[R214] Li, H., & Gnedin, O. Y. 2019, MNRAS, 486, 4030 [NASA ADS] [CrossRef] [Google Scholar]

[R215] Li, H., Gnedin, O. Y., & Gnedin, N. Y. 2018, ApJ, 861, 107 [Google Scholar]

[R216] Li, H., Gnedin, O. Y., Gnedin, N. Y., et al. 2017, ApJ, 834, 69 [NASA ADS] [CrossRef] [Google Scholar]

[R217] Li, H., Vogelberger, M., Bryan, G. L., et al. 2022, MNRAS, 514, 265 [Google Scholar]

[R218] Lim, S., & Lee, M. G. 2015, ApJ, 804, 123 [NASA ADS] [CrossRef] [Google Scholar]

[R219] Lin, C. C., & Shu, F. H. 1966, PNAS, 55, 229 [CrossRef] [Google Scholar]

[R220] Lomelí-Núñez, L., Mayya, Y. D., Rodríguez-Merino, L. H., et al. 2024, MNRAS, 528, 1445 [Google Scholar]

[R221] Mackey, A. D., Wilkinson, M. I., Davies, M. B., & Gilmore, G. F. 2008, MNRAS, 386, 65 [NASA ADS] [CrossRef] [Google Scholar]

[R222] Mai, Y., Croom, S. M., Wisnioski, E., et al. 2024, MNRAS, 533, 15 [Google Scholar]

[R223] Maji, M., Zhu, Q., Li, Y., et al. 2017, ApJ, 844, 108 [NASA ADS] [CrossRef] [Google Scholar]

[R224] Mapelli, M., & Bressan, A. 2013, MNRAS, 430, 3120 [Google Scholar]

[R225] Maraston, C., Bastian, N., Saglia, R. P., et al. 2004, A&A, 416, 467 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[R226] Martin, N. F., Venn, K. A., Aguado, D. S., et al. 2022, Nature, 601, 45 [NASA ADS] [CrossRef] [Google Scholar]

[R227] Maschmann, D., Lee, J. C., Thilker, D. A., et al. 2024, ApJSS, 273, 14 [Google Scholar]

[R228] Matzner, C. D., & McKee, C. F. 2000, ApJ, 545, 364 [NASA ADS] [CrossRef] [Google Scholar]

[R229] McCrady, N., & Graham, J. R. 2007, ApJ, 663, 844 [NASA ADS] [Google Scholar]

[R230] McKee, C. F., & Krumholz, M. R. 2010, ApJ, 709, 308 [NASA ADS] [CrossRef] [Google Scholar]

[R231] McKee, C. F., & Tan, J. C. 2003, ApJ, 585, 850 [Google Scholar]

[R232] Meibom, S., Grundahl, F., Clausen, J. V., et al. 2009, AJ, 137, 5086 [Google Scholar]

[R233] Merritt, D., Piatek, S., Portegies Zwart, S., & Hemsendorf, M. 2004, ApJ, 608, L25 [NASA ADS] [CrossRef] [Google Scholar]

[R234] Metha, B., Trenti, M., & Chu, T. 2021, MNRAS, 508, 489 [NASA ADS] [CrossRef] [Google Scholar]

[R235] Metha, B., Trenti, M., Battisti, A., & Chu, T. 2024, MNRAS, 529, 104 [Google Scholar]

[R236] Meynet, G., Maeder, A., Schaller, G., Schaerer, D., & Charbonnel, C. 1994, A&ASS, 103, 97 [Google Scholar]

[R237] Minniti, D., Geisler, D., Alonso-García, J., et al. 2017, ApJ, 849, L24 [CrossRef] [Google Scholar]

[R238] Monaco, P., Fontanot, F., & Taffoni, G. 2007, MNRAS, 375, 1189 [NASA ADS] [CrossRef] [Google Scholar]

[R239] Mowla, L., Iyer, K. G., Desprez, G., et al. 2022, ApJ, 937, L35 [NASA ADS] [CrossRef] [Google Scholar]

[R240] Mowla, L., Iyer, K., Asada, Y., et al. 2024, Nature, 636, 332 [Google Scholar]

[R241] Muratov, A. L., & Gnedin, O. Y. 2010, ApJ, 718, 1266 [Google Scholar]

[R242] Murphy, G. G., Yates, R. M., & Mohamed, S. S. 2022, MNRAS, 510, 1945 [Google Scholar]

[R243] Myeong, G. C., Evans, N. W., Belokurov, V., Sanders, L. J., & Koposov, S. E. 2018, ApJ, 863, L28 [NASA ADS] [CrossRef] [Google Scholar]

[R244] Myers, P. C., & Goodman, A. A. 1988, ApJ, 329, 392 [Google Scholar]

[R245] Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563 [Google Scholar]

[R246] Newton, O., Davies, J. J., Pfeffer, J., et al. 2025, MNRAS, 542, 591 [Google Scholar]

[R247] Nordlund, Å., & Padoan, P. 1999, in Interstellar Turbulence, eds. J. Franco, & A. Carraminana (Cambridge: Cambridge University Press), 218 [Google Scholar]

[R248] Ohlson, D., Seth, A. C., Gallo, E., Baldassare, V. F., & Greene, J. E. 2023, AJ, 23 [Google Scholar]

[R249] Oklopčić, A., Hopkins, P. F., Feldmann, R., et al. 2017, MNRAS, 465, 952 [CrossRef] [Google Scholar]

[R250] Önehag, A., Gustafsson, B., & Korn, A. 2014, A&A, 562, A102 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[R251] Ono, Y., Ouchi, M., Harikane, Y., et al. 2025, ApJ, 991, 222 [Google Scholar]

[R252] Ostriker, E. C. 1999, ApJ, 513, 252 [Google Scholar]

[R253] Ostriker, E. C., Stone, J. M., & Gammie, C. F. 2001, ApJ, 546, 980 [Google Scholar]

[R254] Pace, A. B. 2025, OJA, 8, 142 [Google Scholar]

[R255] Padoan, P., & Nordlund, Å. 2002, ApJ, 576, 870 [NASA ADS] [CrossRef] [Google Scholar]

[R256] Padoan, P., & Nordlund, Å. 2011, ApJ, 730, 40 [Google Scholar]

[R257] Padoan, P., Jimenez, R., & Jones, B. 1997, MNRAS, 185, 711 [Google Scholar]

[R258] Parente, M., Ragone-Figueroa, C., Granato, G. L., & Lapi, A. 2023, MNRAS, 521, 6105 [NASA ADS] [CrossRef] [Google Scholar]

[R259] Parente, M., Ragone-Figueroa, C., López, P., et al. 2024, ApJ, 966, 154 [NASA ADS] [CrossRef] [Google Scholar]

[R260] Pasquali, A., Bik, A., Zibetti, S., et al. 2011, AJ, 141, 132 [NASA ADS] [CrossRef] [Google Scholar]

[R261] Pastorello, N., Forbes, D. A., Usher, C., et al. 2015, MNRAS, 451, 2625 [Google Scholar]

[R262] Paunzen, E., & Netopil, M. 2006, MNRAS, 371, 1641 [NASA ADS] [CrossRef] [Google Scholar]

[R263] Pechetti, R., Seth, A. C., Neumayer, N., et al. 2020, ApJ, 900, 32 [NASA ADS] [CrossRef] [Google Scholar]

[R264] Peng, E. W., Jordán, A., Côté, P., et al. 2006, ApJ, 639, 95 [Google Scholar]

[R265] Peng, E. W., Jordán, A., Côté, P., et al. 2008, ApJ, 681, 197 [NASA ADS] [CrossRef] [Google Scholar]

[R266] Perrett, K. M., Bridges, T. J., Hanes, D. A., et al. 2002, AJ, 123, 2490 [Google Scholar]

[R267] Pfeffer, J., Kruijssen, J. M. D., Crain, R. A., & Bastian, N. 2018, MNRAS, 475, 4309 [Google Scholar]

[R268] Pfeffer, J., Bastian, N., Crain, R. A., et al. 2019a, MNRAS, 487, 4550 [Google Scholar]

[R269] Pfeffer, J., Bastian, N., Kruijssen, J. M. D., et al. 2019b, MNRAS, 490, 1714 [Google Scholar]

[R270] Pfeffer, J., Kruijssen, J. M. D., Bastian, N., Crain, R. A., & Trujillo-Gomez, S. 2023, MNRAS, 519, 5384 [Google Scholar]

[R271] Pfeffer, J., Janssens, S. R., Buzzo, M. L., et al. 2024, MNRAS, 429, 4914 [Google Scholar]

[R272] Pfeffer, J., Forbes, D. A., Romanowsky, A. J., et al. 2025, MNRAS, 536, 1878 [Google Scholar]

[R273] Planck Collaboration XVI. 2014, A&A, 571, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[R274] Plummer, H. C. 1911, MNRAS, 71, 460 [Google Scholar]

[R275] Poggio, E., Recio-Blano, A., Palicio, P. A., et al. 2022, A&A, 666, L4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[R276] Polkas, M., Bonoli, S., Bortolas, E., et al. 2024, A&A, 689, A204 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[R277] Portegies Zwart, S. F., McMillan, S. L. W., & Gieles, M. 2010, ARAA, 48, 431 [Google Scholar]

[R278] Prieto, J. L., & Gnedin, O. Y. 2008, ApJ, 689, 919 [Google Scholar]

[R279] Rafelski, M., & Zaritsky, D. 2005, AJ, 129, 2701 [NASA ADS] [CrossRef] [Google Scholar]

[R280] Rafikov, R. R. 2001, MNRAS, 323, 445 [NASA ADS] [CrossRef] [Google Scholar]

[R281] Ramírez-Alegría, S., Borissova, J., Chené, A. B., et al. 2014, A&A, 564, L9 [Google Scholar]

[R282] Randriamanakoto, Z., Väisänen, P., Ryder, S. D., & Ranaivomanana, P. 2019, MNRAS, 482, 2530 [Google Scholar]

[R283] Reina-Campos, M., & Kruijssen, J. M. D. 2017, MNRAS, 469, 1282 [CrossRef] [Google Scholar]

[R284] Reina-Campos, M., Kruijssen, J. M. D., Pfeffer, J. L., Bastian, N., & Crain, R. A. 2019, MNRAS, 486, 5838 [NASA ADS] [CrossRef] [Google Scholar]

[R285] Reina-Campos, M., Hughes, M. E., Kruijssen, J. M. D., et al. 2020, MNRAS, 493, 3422 [Google Scholar]

[R286] Reina-Campos, M., Keller, B. W., Kruijssen, J. M. D., et al. 2022a, MNRAS, 517, 3144 [NASA ADS] [CrossRef] [Google Scholar]

[R287] Reina-Campos, M., Trujillo-Gomez, S., Deason, A. J., et al. 2022b, MNRAS, 513, 3925 [NASA ADS] [CrossRef] [Google Scholar]

[R288] Reina-Campos, M., Sills, A., & Bichon, G. 2023a, MNRAS, 524, 968 [NASA ADS] [CrossRef] [Google Scholar]

[R289] Reina-Campos, M., Trujillo-Gomez, S., Pfeffer, J. L., et al. 2023b, MNRAS, 521, 6368 [NASA ADS] [CrossRef] [Google Scholar]

[R290] Reina-Campos, M., Gnedin, O. Y., Sills, A., & Li, H. 2025, ApJ, 978, 15 [Google Scholar]

[R291] Renaud, F. 2018, New Astron. Rev., 81, 1 [CrossRef] [Google Scholar]

[R292] Renaud, F. 2020, IAU Symp., 351, 40 [Google Scholar]

[R293] Renaud, F., Gieles, M., & Boily, C. M. 2011, MNRAS, 418, 759 [NASA ADS] [CrossRef] [Google Scholar]

[R294] Rodruck, M., Charlton, J., Borthakur, S., et al. 2023, MNRAS, 526, 2341 [Google Scholar]

[R295] Romeo, A. B., & Wiegert, J. 2011, MNRAS, 416, 1191 [NASA ADS] [CrossRef] [Google Scholar]

[R296] Rosolowsky, E., & Blitz, L. 2005, ApJ, 623, 826 [NASA ADS] [CrossRef] [Google Scholar]

[R297] Rostami-Shirazi, A., Zonoozi, A. H., Haghi, H., & Rabiee, M. 2024, MNRAS, 535, 3489 [Google Scholar]

[R298] Ryon, J. E., Adamo, A., Bastian, N., et al. 2014, AJ, 148, 33 [Google Scholar]

[R299] Ryon, J. E., Bastian, N., Adamo, A., et al. 2015, MNRAS, 452, 525 [CrossRef] [Google Scholar]

[R300] Ryon, J. E., Gallagher, J. S., Smith, L. J., et al. 2017, ApJ, 841, 92 [Google Scholar]

[R301] Safronov, V. S. 1960, Annales d’Astrophysique, 23, 979 [NASA ADS] [Google Scholar]

[R302] Saha, K., tseng, Y.-H., & Taam, R. E. 2010, ApJ, 721, 1878 [Google Scholar]

[R303] Salmon, B., Coe, D., Bradley, L., et al. 2018, ApJ, 864, L22 [Google Scholar]

[R304] Sánchez-Menguiano, L., Sánchez, S. F., Kawata, D., et al. 2016, ApJ, 830, L40 [CrossRef] [Google Scholar]

[R305] Sánchez-Menguiano, L., Sánchez, S. F., Pérez, I., et al. 2018, A&A, 609, A119 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[R306] Sanders, R. L., Shapley, A. E., Jones, T., et al. 2021, ApJ, 914, 19 [CrossRef] [Google Scholar]

[R307] Sattari, Z., Mobasher, B., Chartab, N., et al. 2023, ApJ, 951, 147 [NASA ADS] [CrossRef] [Google Scholar]

[R308] Schaller, G., Schaerer, D., Meynet, G., & Maeder, A. 1992, A&ASS, 96, 269 [Google Scholar]

[R309] Schaye, J., Crain, R. A., Bower, R. G., et al. 2015, MNRAS, 446, 521 [Google Scholar]

[R310] Schuster, K. F., Kramer, C., Hitschfeld, M., Garcia-Burillo, S., & Mookerjea, B. 2007, A&A, 461, 143 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[R311] Sesto, L. A., Faifer, F. R., Smith Castelli, A. V., Forte, J. C., & Escudero, C. G. 2018, MNRAS, 479, 478 [NASA ADS] [Google Scholar]

[R312] Shibuya, T., Ouchi, M., Kubo, M., & Harikane, Y. 2016, ApJ, 821, 72 [NASA ADS] [CrossRef] [Google Scholar]

[R313] Silva-Villa, E., & Larsen, S. S. 2011, A&A, 529, A25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[R314] Somerville, R. S., & Davé, R. 2015, ARAA, 53, 51 [Google Scholar]

[R315] Somerville, R. S., Hopkins, P. F., Cox, T. J., Robertson, B. E., & Hernquist, L. 2008, MNRAS, 391, 481 [NASA ADS] [CrossRef] [Google Scholar]

[R316] Spengler, C., Côté, P., Roediger, J., et al. 2017, ApJ, 849, 55 [NASA ADS] [CrossRef] [Google Scholar]

[R317] Spina, L., Magrini, L., & Cunha, K. 2022, Universe, 8, 87 [NASA ADS] [CrossRef] [Google Scholar]

[R318] Spinoso, D., Bonoli, S., Valiante, R., Schneider, R., & Izquierdo-Villalba, D. 2023, MNRAS, 518, 4672 [Google Scholar]

[R319] Spitler, L. R., & Forbes, D. A. 2009, MNRAS, 392, L1 [Google Scholar]

[R320] Spitzer, L. 1940, MNRAS, 100, 396 [Google Scholar]

[R321] Spitzer, L. 1987, Dynamical Evolution of Globular Clusters (Princeton: Princeton University Press) [Google Scholar]

[R322] Springel, V., White, S. D. M., Tormen, G., & Kauffmann, G. 2001, MNRAS, 328, 726 [Google Scholar]

[R323] Springel, V., White, S. D. M., Jenkins, A., et al. 2005, Nature, 435, 629 [Google Scholar]

[R324] Strader, J., Brodie, J. P., Cenarro, A. J., Beasley, M. A., & Forbes, D. A. 2005, AJ, 130, 1315 [NASA ADS] [CrossRef] [Google Scholar]

[R325] Straižys, V., Maskoliūnas, M., Boyle, R. P., et al. 2014, MNRAS, 437, 1628 [Google Scholar]

[R326] Sun, J., Leroy, A. K., Rosolowsky, E., et al. 2022, AJ, 164, 43 [NASA ADS] [CrossRef] [Google Scholar]

[R327] Tanaka, T., & Haiman, Z. 2009, ApJ, 696, 1798 [NASA ADS] [CrossRef] [Google Scholar]

[R328] Thilker, D. A., Lee, J. C., Whitmore, B. C., et al. 2025, ApJSS, 280, 1 [Google Scholar]

[R329] Tonini, C. 2013, ApJ, 762, 39 [NASA ADS] [CrossRef] [Google Scholar]

[R330] Toomre, A. 1964, ApJ, 139, 1217 [Google Scholar]

[R331] Torrey, P., Vogelsberger, M., Marinacci, F., et al. 2019, MNRAS, 484, 5587 [NASA ADS] [Google Scholar]

[R332] Tremonti, C. A., Heckman, T. M., Kauffmann, G., et al. 2004, ApJ, 613, 898 [Google Scholar]

[R333] Trujillo-Gomez, S., Kruijssen, J. M. D., Reina-Campos, M., et al. 2021, MNRAS, 503, 31 [NASA ADS] [CrossRef] [Google Scholar]

[R334] Untzaga, J., Bonoli, S., Izquierdo-Villalba, D., Mezcua, M., & Spinoso, D. 2024, MNRAS, 535, 3293 [Google Scholar]

[R335] Usher, C., Forbes, D. A., Brodie, J. P., et al. 2012, MNRAS, 426, 1475 [Google Scholar]

[R336] Usher, C., Forbes, D. A., Spitler, L. R., et al. 2013, MNRAS, 436, 1172 [Google Scholar]

[R337] Usher, C., Caldwell, N., & Cabrera-Ziri, I. 2024, MNRAS, 528, 6010 [NASA ADS] [CrossRef] [Google Scholar]

[R338] Valenzuela, L. M., Moster, B. P., Remus, R.-S., O’Leary, J. A., & Burkert, A. 2021, MNRAS, 505, 5815 [NASA ADS] [CrossRef] [Google Scholar]

[R339] Valenzuela, L. M., Forbes, D. A., & Remus, R.-S. 2025, MNRAS, 537, 306 [CrossRef] [Google Scholar]

[R340] van der Kruit, P. C., & Freeman, K. C. 2011, ARAA, 49, 301 [Google Scholar]

[R341] Vani, A., Ayromlou, M., Kauffmann, G., & Springel, V. 2025, MNRAS, 536, 777 [NASA ADS] [Google Scholar]

[R342] Vanzella, E., Calura, F., Meneghetti, M., et al. 2017, MNRAS, 467, 4303 [Google Scholar]

[R343] Vanzella, E., Calura, F., Meneghetti, M., et al. 2019, MNRAS, 483, 3618 [Google Scholar]

[R344] Vanzella, E., Claeyssens, A., Welch, B., et al. 2023, ApJ, 945, 53 [NASA ADS] [CrossRef] [Google Scholar]

[R345] Vazquez-Semadeni, E. 1994, ApJ, 423, 681 [Google Scholar]

[R346] Vázquez-Semadeni, J. S. E., Chappell, D., & Passot, T. 1998, ApJ, 504, 835 [CrossRef] [Google Scholar]

[R347] Veljanoski, J., Ferguson, A. M. N., Mackey, A. D., et al. 2015, MNRAS, 452, 320 [NASA ADS] [CrossRef] [Google Scholar]

[R348] Vijayan, A. P., Scott, S. J., Thomas, P. A., et al. 2019, MNRAS, 489, 4072 [Google Scholar]

[R349] Villaume, A., Romanowsky, A. J., Brodie, J., & Strader, J. 2019, ApJ, 879, 45 [NASA ADS] [CrossRef] [Google Scholar]

[R350] Wang, K., & Peng, Y. 2024, A&A, 12 [Google Scholar]

[R351] Wang, E., Wang, H., Mo, H., et al. 2018, ApJ, 864, 51 [NASA ADS] [CrossRef] [Google Scholar]

[R352] Watson, G. N. 1944, A Treatise on the Theory of Bessel Functions, The Mathematical Gazette (Cambridge: Cambridge University Press) [Google Scholar]

[R353] West, M. J., Côté, P., Jones, C., Forman, W., & Marzke, R. O. 1995, ApJ, 453, L77 [NASA ADS] [CrossRef] [Google Scholar]

[R354] West, M. J., Côté, P., Marzke, R. O., & Jordán, A. 2004, Nature, 427, 31 [Google Scholar]

[R355] Westfall, K. B., Andersen, D. R., Bershady, M. A., et al. 2014, ApJ, 785, 43 [Google Scholar]

[R356] White, S. D. M., & Frenk, C. S. 1991, ApJ, 379, 52 [Google Scholar]

[R357] White, S. D. M., & Rees, M. J. 1978, MNRAS, 183, 341 [Google Scholar]

[R358] Whitmore, B. C., Chandar, R., Bowers, A. S., et al. 2014, AJ, 147, 78 [NASA ADS] [CrossRef] [Google Scholar]

[R359] Yates, R. M., Henriques, B. M. B., Thomas, P. A., et al. 2013, MNRAS, 435, 3500 [Google Scholar]

[R360] Yates, R. M., Thomas, P. A., & Henriques, B. M. B. 2017, MNRAS, 464, 3169 [NASA ADS] [CrossRef] [Google Scholar]

[R361] Yates, R. M., Henriques, B. M. B., Fu, J., et al. 2021, MNRAS, 503, 4474 [NASA ADS] [CrossRef] [Google Scholar]

[R362] Yates, R. M., Hendriks, D., Vijayan, A. P., et al. 2024, MNRAS, 527, 6292 [Google Scholar]

[R363] Zaritsky, D. 2022, MNRAS, 513, 2609 [NASA ADS] [CrossRef] [Google Scholar]

[R364] Zentner, A. R., & Bullock, J. S. 2003, ApJ, 598, 49 [NASA ADS] [CrossRef] [Google Scholar]

[R365] Zhang, Q., & Fall, S. M. 1999, ApJ, 527, L81 [CrossRef] [Google Scholar]

[R366] Zhong, W., Fu, J., Sharma, P., Shen, S., & Yates, R. M. 2023a, MNRAS, 519, 4344 [Google Scholar]

[R367] Zhong, W., Fu, J., Shen, S., & Yuan, F. 2023b, RAA, 23, 075004 [Google Scholar]

[R368] Zhou, L., Federrath, C., Yuan, T., et al. 2017, MNRAS, 470, 4573 [Google Scholar]

Massive star clusters in the semi-analytical galaxy formation model L-Galaxies 2020

1 Introduction

2 Model description

2.1 The L-Galaxies semi-analytical galaxy formation model

2.2 Gravitational potential

2.3 Star cluster formation

2.3.1 Epicyclic frequency

2.3.2 Toomre stability parameter

2.3.3 Bound fraction of star formation

2.4 Initial star cluster properties

2.4.1 Stellar mass

2.4.2 Initial galactocentric distances

2.4.3 Half-mass radius

2.4.4 Tidal radius

2.4.5 Metallicity

2.5 Star cluster evolution

2.5.1 Mass loss

2.5.2 Radial expansion

2.5.3 Galactocentric migration

2.5.4 Re-distribution during galaxy mergers

3 Results

3.1 MV-SFR relationship

3.2 Stellar masses

3.3 Half-mass radii

3.4 Galactocentric distances

3.5 Metallicity distributions

3.5.1 Mean cluster metallicity

3.5.2 Bimodality

3.6 Caveats

3.6.1 Axisymmetric structures of galaxies

3.6.2 Impact of stellar-mass black holes

3.6.3 Galaxy-galaxy interactions

4 Conclusions

Data availability

Acknowledgements

References

Appendix A Model variations

Appendix A.1 Velocity dispersion of the cold gas

Appendix A.2 Initial half-mass radius

Appendix A.3 Cluster initial mass function

All Tables

All Figures

3.1 M_V-SFR relationship