Globular cluster metallicity distributions in the E-MOSAICS simulations Open Access

ABSTRACT

The metallicity distributions of globular cluster (GC) systems in galaxies are a critical test of any GC formation scenario. In this work, we investigate the predicted GC metallicity distributions of galaxies in the MOdelling Star cluster population Assembly In Cosmological Simulations within EAGLE (E-MOSAICS) simulation of a representative cosmological volume (L = 34.4 comoving Mpc). We find that the predicted GC metallicity distributions and median metallicities from the fiducial E-MOSAICS GC formation model agree well the observed distributions, except for galaxies with masses |$M_\ast \sim 2 \times 10^{10} \, \rm {M}_{\odot }$|⁠, which contain an overabundance of metal-rich GCs. The predicted fraction of galaxies with bimodal GC metallicity distributions (37 ± 2 per cent in total; 45 ± 7 per cent for |$M_\ast \gt 10^{10.5} \, \rm {M}_{\odot }$|⁠) is in good agreement with observed fractions (⁠|$44^{+10}_{-9}$| per cent), as are the mean metallicities of the metal-poor and metal-rich peaks. We show that, for massive galaxies (⁠|$M_\ast \gt 10^{10} \, \rm {M}_{\odot }$|⁠), bimodal GC distributions primarily occur as a result of cluster disruption from initially-unimodal distributions, rather than as a result of cluster formation processes. Based on the distribution of field stars with GC-like abundances in the Milky Way, we suggest that the bimodal GC metallicity distribution of Milky Way GCs also occurred as a result of cluster disruption, rather than formation processes. We conclude that separate formation processes are not required to explain metal-poor and metal-rich GCs, and that GCs can be considered as the surviving analogues of young massive star clusters that are readily observed to form in the local Universe today.

1 INTRODUCTION

Globular clusters (GCs) are one of the most common types of stellar system in the Universe, with most galaxies with stellar masses larger than |$10^9 \, \rm {M}_{\odot }$| hosting GC populations (Harris 1991; Brodie & Strader 2006). However, their formation mechanism is still under debate (see Kruijssen 2014; Forbes et al. 2018, for recent reviews), with proposed scenarios broadly falling into two classes: special conditions for GC formation existed in the Early Universe (e.g. Peebles & Dicke 1968; Peebles 1984; Fall & Rees 1985; Rosenblatt, Faber & Blumenthal 1988; Trenti, Padoan & Jimenez 2015; Mandelker et al. 2018; Madau et al. 2020); or, a universal formation mechanism explains both GCs and young star clusters (e.g. Ashman & Zepf 1992; Harris & Pudritz 1994; Kravtsov & Gnedin 2005; Kruijssen 2015; Li et al. 2017; Pfeffer et al. 2018; Ma et al. 2020; Reina-Campos et al. 2022b; Grudić et al. 2023).

One of the most peculiar properties of the GC systems of massive galaxies is their apparent colour and/or metallicity bimodality. Since the discovery of metallicity bimodality in the Milky Way’s GC population (Harris & Canterna 1979; Freeman & Norris 1981; Zinn 1985), it has been a major focus in studies of extragalactic GC systems (for reviews, see Ashman & Zepf 1998; Harris 2001; Brodie & Strader 2006). Colour bimodality (or multimodality) appears commonly in massive galaxies (e.g. Gebhardt & Kissler-Patig 1999; Larsen et al. 2001), with a contribution of red GCs that decreases towards fainter galaxies (Peng et al. 2006), though the interpretation as two distinct populations becomes complicated where internal colour dispersions (inferred through Gaussian fits) become large (Harris et al. 2017).¹ The two populations also often (e.g. Zinn 1985; Arnold et al. 2011), though not always (e.g. Richtler et al. 2004; Dolfi et al. 2021), show evidence of distinct kinematic signatures, providing further evidence for being distinct populations.

Due to non-linear colour–metallicity relations of GC systems, colour bimodality does not necessarily imply metallicity bimodality (Peng et al. 2006; Yoon, Yi & Lee 2006; Conroy, Gunn & White 2009; Usher et al. 2012). The GC colour–metallicity relations may also vary from galaxy to galaxy (Usher et al. 2015; Sesto et al. 2018; Villaume et al. 2019), possibly caused by age or abundance differences, making the inference of metallicity distributions from colour distributions less than straightforward. Therefore, spectroscopic studies are necessary to confirm the intrinsic metallicity distributions of GC populations.

The more expensive nature of spectroscopy, compared to photometry, means that fewer systematic spectroscopic studies have been undertaken on the metallicity distributions of the GC populations of galaxies. This is especially the case for typical star-forming galaxies, since most studies focus on quiescent early-type galaxies, where extinction and reddening from dust are less of an issue. While NGC 3115 shows one of the best examples of a bimodal GC metallicity distribution (Brodie et al. 2012), other galaxies often show more diverse distributions (Beasley et al. 2008; Usher et al. 2012; Fahrion et al. 2020b). In the Local Group, the GC metallicity distribution of M31 appears to be broad and unimodal, in contrast with the bimodal distribution in the Milky Way (Caldwell et al. 2011). Thus, though the bimodality of GC colour distributions may be common, this is not necessarily a universal feature of GC metallicity distributions. Any complete scenario for GC formation should therefore explain the origin of the diversity of GC metallicity distributions.

Many scenarios for GC formation have been developed and interpreted through the lens of bimodality. Scenarios often invoke separate formation mechanisms and/or truncated formation epochs in order to explain metal-poor and metal-rich GC populations (e.g. Zepf & Ashman 1993; Forbes, Brodie & Grillmair 1997; Beasley et al. 2002; Griffen et al. 2010). However a single, continuous formation mechanism may also produce GC bimodality (Kravtsov & Gnedin 2005; Bekki et al. 2008; Muratov & Gnedin 2010; Li & Gnedin 2014, 2019; Kruijssen 2015; Choksi, Gnedin & Li 2018; Keller et al. 2020). Regardless of the actual formation mechanism of GCs, it is thought to be a natural outcome of hierarchical galaxy formation (at least in massive galaxies) that the metal-rich and metal-poor GCs largely represent in situ and ex situ formation within their host galaxy, respectively (Côté, Marzke & West 1998; Hilker, Infante & Richtler 1999; Tonini 2013; Katz & Ricotti 2014; Renaud, Agertz & Gieles 2017; Forbes & Remus 2018; El-Badry et al. 2019; Kruijssen et al. 2019a,b). In this way, the metallicity distributions of GC systems may be intimately tied to the formation and assembly history of their host galaxy, with differing assembly histories resulting in differing GC metallicity distributions (e.g. Tonini 2013; Kruijssen et al. 2019a). Often neglected in studies of GC metallicity distributions, cluster mass loss and disruption may also significantly alter the distributions. For instance, Kruijssen (2015) showed that mass-loss mechanisms are expected to be more effective for higher metallicity GCs. Therefore, present-day distributions could be significantly different from the initial distributions, complicating the inference of formation scenarios from present-day GC population properties.

Closely related to GC metallicity distributions are the specific frequencies of GCs (the number of GCs relative to the total luminosity, mass or number of field stars) as a function of metallicity in galaxies. Observations show that the specific frequency is a strong function of metallicity and decreases with increasing metallicity; i.e. more stars are found in GCs at low metallicities (⁠|$[\rm {Fe}/\rm {H} ]\sim -2$|⁠) than at higher metallicities (Durrell, Harris & Pritchet 2001; Harris & Harris 2002; Beasley et al. 2008; Gratton, Carretta & Bragaglia 2012; Kruijssen 2015; Lamers et al. 2017). This relation could be imparted either at formation (more efficient formation of low metallicity clusters), by cluster disruption (higher mass-loss rates for higher metallicity clusters), or some combination of both, and presents an important diagnostic for models of GC system formation and evolution (see Lamers et al. 2017, for a comprehensive discussion).

In this work, we investigate the metallicity distributions and metallicity-dependent specific frequencies of GCs in the MOdelling Star cluster population Assembly In Cosmological Simulations within EAGLE (E-MOSAICS; Pfeffer et al. 2018; Kruijssen et al. 2019a) simulation of galaxy and star cluster formation within a periodic cosmological volume of side length L = 34.4 comoving Mpc. E-MOSAICS couples models for star cluster formation and evolution to the Evolution and Assembly of GaLaxies and their Environments (EAGLE) galaxy formation model (Schaye et al. 2015; Crain et al. 2015), under the assumption that both young and old star clusters form and evolve following the same physical mechanisms. We aim to test whether a common formation mechanism for young and old clusters can explain the diversity of metallicity distributions found for GC systems of observed galaxies.

This paper is organised as follows. In Section 2 we briefly describe the E-MOSAICS model and our GC selection methods for the simulation. Section 3 presents the results from this study, including predictions for GC metallicity distributions, specific frequency–metallicity relations and dependence on the GC formation model. In Section 4 we discuss the origin of bimodal GC metallicity distributions and connection to GC formation scenarios. Finally, in Section 5 we summarise the conclusions of this work.

2 METHODS

In this section we briefly describe the E-MOSAICS and EAGLE models, the periodic volume simulation analysed in this work and our selection of GCs from the simulation.

2.1 E-MOSAICS simulations

E-MOSAICS (Pfeffer et al. 2018; Kruijssen et al. 2019a) is a suite of cosmological, hydrodynamical simulations that couple the MOSAICS (Kruijssen et al. 2011; Pfeffer et al. 2018) model for star cluster formation and evolution to the EAGLE galaxy formation model (Crain et al. 2015; Schaye et al. 2015). The simulations were performed with a highly-modified version of the N-body, smooth particle hydrodynamics code gadget 3 (Springel 2005) and the EAGLE model includes subgrid routines for radiative cooling (Wiersma, Schaye & Smith 2009a), star formation (Schaye & Dalla Vecchia 2008), stellar evolution (Wiersma et al. 2009b), the seeding and growth of supermassive black holes (Rosas-Guevara et al. 2015), and feedback from star formation and black hole growth (Booth & Schaye 2009). In EAGLE, the feedback efficiencies for stellar and black hole feedback are calibrated such that simulations of representative galaxy populations reproduce the galaxy stellar mass function, galaxy sizes and black hole masses at z ≈ 0 (Crain et al. 2015). Simulations using the EAGLE model have been well studied and shown to broadly reproduce many features of the evolving galaxy population, including the evolution of the galaxy stellar mass function (Furlong et al. 2015) and galaxy sizes (Furlong et al. 2017), galaxy star formation rates and colours (Furlong et al. 2015; Trayford et al. 2017), cold gas properties (Lagos et al. 2015; Crain et al. 2017), galaxy morphologies (Bignone et al. 2020), and (particularly relevant for this work) the galaxy mass-metallicity relation (Schaye et al. 2015, with the high-resolution ‘recalibrated’ model used in this work).

In the MOSAICS model, star clusters are treated as subgrid components of stellar particles, such that they adopt the properties of their host particle (positions, velocities, ages, abundances). Star cluster formation is described by two functions: the fraction of stars formed in bound clusters (i.e. the cluster formation efficiency, CFE, Bastian 2008) and the shape of the initial cluster mass function (a power law or Schechter 1976 function, i.e. power law with a high mass exponential truncation |$M_{\rm {c,\ast}}$|⁠). For both the power law and Schechter initial mass functions the initial power law index is set to be −2, based on observations of young clusters (e.g. Zhang & Fall 1999; Bik et al. 2003; Gieles et al. 2006b; McCrady & Graham 2007). In each case, the CFE or |$M_{\rm c,\ast }$| may be constant or vary with the properties of the local environment at the time of star formation (see below). Each newly-formed star particle may (stochastically) form a fraction of its mass in clusters (i.e. the CFE times the particle mass) and clusters are sampled stochastically from the initial mass function (such that the subgrid clusters may be more massive than the stellar mass of the host particle). Thus the total cluster and field star mass is only conserved for an ensemble of star particles (for details of the method, see Pfeffer et al. 2018). Following formation, star clusters may lose mass by stellar evolution (following the EAGLE model), two-body relaxation depending on the local tidal field strength (Lamers et al. 2005; Kruijssen et al. 2011; with an additional constant term to account for isolated clusters, following Gieles & Baumgardt 2008) and tidal shocks from rapidly changing tidal fields (Gnedin, Hernquist & Ostriker 1999; Prieto & Gnedin 2008; Kruijssen et al. 2011). Dynamical friction is treated in post-processing at every snapshot and assumed to completely remove clusters when the dynamical friction time-scale is less than the cluster age (assuming they merge to the centre of their host galaxy, see Pfeffer et al. 2018).

In the E-MOSAICS suite we consider four main variations of the cluster formation model, i.e. two variations each (constant or environmentally varying) for the CFE and initial cluster mass function. Comparison of the four models enables the most critical aspect of cluster formation to be determined for a particular observable. The environmentally-varying CFE is determined by the Kruijssen (2012) model, which scales as a function of the natal gas pressure such that higher pressures result in higher CFE (see fig. 3 in Pfeffer et al. 2018). The environmentally-varying mass function has an exponential truncation (⁠|$M_{\rm c,\ast }$|⁠) which varies with local gas and dynamical properties according to the model of Reina-Campos & Kruijssen (2017), such that |$M_{\rm c,\ast }$| increases with local gas pressure, except where limited by high Coriolis or centrifugal forces (i.e. near the centres of galaxies). Thus the four cluster formation models are:

• Fiducial: Both the CFE and |$M_{\rm c,\ast }$| are environmentally dependent. The fiducial model reproduces observed scaling relations of young star clusters (Pfeffer et al. 2019).

• CFE only: Environmentally-varying CFE and pure power-law mass function (i.e. |$M_{\rm c,\ast }= \infty$|⁠).

• |$M_{\rm c,\ast }$| only: Constant CFE = 0.1 and environmentally-varying |$M_{\rm c,\ast }$|⁠.

• Constant formation: Constant CFE = 0.1 and power-law mass function (⁠|$M_{\rm c,\ast }= \infty$|⁠).

The fiducial E-MOSAICS model has previously been shown to produce GC populations consistent with many observed relations, including the fraction of stars contained in GCs (Bastian et al. 2020), GC system radial distributions (Reina-Campos et al. 2022a), the high-mass truncation of GC mass functions (Hughes et al. 2022), the ‘blue tilt’ GC colour distributions (Usher et al. 2018), and the age–metallicity relations of GC systems (Kruijssen et al. 2019b, 2020; Horta et al. 2021a). The alternative formation models (CFE only, |$M_{\rm c,\ast }$| only, constant formation) generally fail to simultaneously reproduce particular aspects of star cluster populations (Usher et al. 2018; Pfeffer et al. 2019; Reina-Campos et al. 2019; Bastian et al. 2020; Hughes et al. 2022). As discussed in Pfeffer et al. (2018) and Kruijssen et al. (2019a), in Milky Way-mass haloes the E-MOSAICS simulations over-predict the number of low-mass and high-metallicity GCs. This issue is postulated to be a result of a lack of a dense substructure in the interstellar medium (ISM) which may disrupt clusters through tidal shocks, as the EAGLE model does not simulate the cold, dense ISM phase. This issue may be remedied in simulations which include models for the cold ISM (e.g. Reina-Campos et al. 2022b).

In this work, we analyse the E-MOSAICS simulation of a periodic volume with side length 34.4 comoving Mpc (first reported in Bastian et al. 2020). In the simulation, all four cluster formation models were performed in parallel, which is possible as the EAGLE galaxy model is independent of the MOSAICS cluster model. The simulation adopts the ‘Recalibrated’ EAGLE model parameters (Schaye et al. 2015) and initially has 1034³ dark matter and gas particles with masses of |$m_\mathrm{dm} = 1.21 \times 10^6 \, \rm {M}_{\odot }$| and |$m_\mathrm{b} = 2.26 \times 10^5 \, \rm {M}_{\odot }$|⁠, respectively. The simulation adopts parameters consistent with a Planck Collaboration XVI (2014) cosmology (Ω_m = 0.307, |$\Omega _\Lambda = 0.693$|⁠, Ω_b = 0.04825, h = 0.6777, σ₈ = 0.8288). Galaxies (subhaloes) were identified in simulation snapshots by first detecting dark matter structures with the friends-of-friends (FoF) algorithm (Davis et al. 1985) and then identifying bound subhaloes within each FoF group using the subfind (Springel et al. 2001; Dolag et al. 2009) algorithm. The subhalo containing the particle with the minimum gravitational potential is defined as the central galaxy in each FoF group.

2.2 Globular cluster selection

Given the large variations in GC selection between different observational studies/individual galaxies (e.g. due to differing distances of galaxies, background field star surface brightnesses, spectroscopic limitations, etc.), it is difficult to perform systematic comparisons with consistent GC selection. We base our GC selection on the GC luminosity limits of the ACS Virgo/Fornax surveys. In particular, we aim to compare predictions from the simulations against observations of the Virgo cluster from Peng et al. (2006) and the Fornax cluster from Fahrion et al. (2020a, b). Additionally, we select GCs with ages |${\gt} 2 \, \rm {Gyr}$| (i.e. we exclude young star clusters which may undergo significant mass loss through stellar evolution) and metallicities |$[\rm {Fe}/\rm {H} ]\gt -3$| (due to resolution limits, i.e. baryonic particle masses |${\approx}2 \times 10^5 \, \rm {M}_{\odot }$|⁠, very low metallicity particles are poorly sampled in the simulations).

Following table 1 in Jordán et al. (2007, with galaxy stellar masses from Peng et al. 2006), GCs in |$10^9 \, \rm {M}_{\odot }$| galaxies in Virgo are ∼90 per cent complete at the lower apparent magnitude limit of m_g ≈ 26.4 mag. Virgo galaxies with stellar masses |${\approx}10^{11} \, \rm {M}_{\odot }$| have 90 per cent GC completeness at m_g ≈ 24.3 mag. The difference in distance between the Virgo and Fornax clusters implies mass limits around 40 per cent larger in the Fornax cluster (assuming distance moduli of 31.09 and 31.51, respectively; Peng et al. 2006; Blakeslee et al. 2009). Thus we adopt lower GC luminosity limits in the ACS g band with a floor at M_g = −5 and a ceiling at M_g = −7, that scales linearly with the logarithmic stellar mass from M_g = −5 at log (M_*/M_⊙) = 9 to M_g = −7 at log (M_*/M_⊙) = 11. We compare our adopted limits with the Virgo galaxy completeness limits in Fig. 1.

Colours were generated for the simulated GCs using the Flexible Stellar Population Synthesis (fsps) models (Conroy et al. 2009; Conroy & Gunn 2010), assuming simple stellar populations and using the MILES spectral library (Sánchez-Blázquez et al. 2006), Padova isochrones (Girardi et al. 2000; Marigo & Girardi 2007; Marigo et al. 2008), a Chabrier (2003) initial stellar mass function (as assumed in the EAGLE model), and assuming the default fsps parameters. Mass-to-light ratios for the GCs were calculated by linearly interpolating the from the grid in ages and total metallicities (log Z/Z_⊙).

3 RESULTS

3.1 Average metallicity distributions

In Fig. 2 we show the average GC metallicity distributions (i.e. stacked distributions normalized by number of galaxies) predicted in the E-MOSAICS simulations as a function of galaxy mass. For the fiducial comparisons we consider GCs with ages |${\gt}8 \, \rm {Gyr}$|⁠. This age selection is similar to the youngest Milky Way GCs (e.g. Forbes & Bridges 2010) and most similar to the old ages of GCs in cluster galaxies (which undergo early star formation quenching, e.g. Gallazzi et al. 2021). We also show predictions for ages |${\gt}2 \, \rm {Gyr}$| as black dash–dotted lines for reference. The effect of a younger age limit (⁠|$2 \, \rm {Gyr}$|⁠) shows the most difference in galaxies with masses |$8.5 \lesssim \log (M_\ast / \rm {M}_{\odot }) \lt 10$|⁠, which is due to galaxy downsizing (i.e. lower-mass galaxies form more mass at later times, Bower, Lucey & Ellis 1992; Gallazzi et al. 2005).

In the figure we also compare the simulations against observations from the Fornax cluster (Fahrion et al. 2020a, b), SLUGGS survey (Usher et al. 2012, 2015, 2019), and local field galaxies (Milky Way, Harris 1996, 2010; M31, Caldwell et al. 2011; M81, Nantais & Huchra 2010; Large and Small Magellanic Cloud, Horta et al. 2021a; Fornax dwarf spheroidal, Larsen, Brodie & Strader 2012). Where possible, for Local Group galaxies we also apply the same |${\gt}8 \, \rm {Gyr}$| cluster age limit as the simulation (i.e. for the Large and Small Magellanic Clouds). We note that the SLUGGS GC samples are not complete (particularly for the most massive galaxies), and thus only the relative distributions can be compared in these cases. Additionally, we have not attempted to match the environments or morphologies of the observed galaxies, which vary significantly (i.e. from field galaxies to galaxy clusters), and simply compare all galaxies contained in the periodic volume.

In general, the predictions from simulations show a good match to observed GC metallicity distributions. In most cases, the observed distributions for each galaxy mass are fully contained within the galaxy-to-galaxy scatter of the simulations (grey shaded regions), except for the highest mass galaxies (⁠|$M_\ast \gt 10^{11} \, \rm {M}_{\odot }$|⁠) and very low metallicity GCs (⁠|$[\rm {Fe}/\rm {H} ]\lt -2.5$|⁠) which we discuss below. In the 9.5 < log₁₀(M_*/M_⊙) < 10 mass range GCs with metallicities |$[\rm {Fe}/\rm {H} ]\sim -1.5$| appear under-abundant in simulated galaxies by a factor ≈2. Potentially, this could be due to inefficient formation of GCs in such galaxies, which we discuss further in Section 3.3, or simply that our GC selection (Section 2.2) is too strict in comparison to these Fornax cluster galaxies. For Milky Way-mass galaxies [10.5 < log₁₀(M_*/M_⊙) < 11], the simulations predict a factor ≈2.5 times more high-metallicity GCs (⁠|$[\rm {Fe}/\rm {H} ]\sim -0.5$|⁠) than typically observed. This was previously discussed by Kruijssen et al. (2019a) and suggested to be caused by insufficient GC disruption in galactic discs in the simulations (see also Section 2.1).

In the figure we also compare the ‘initial’ GC metallicity distributions (black dotted lines; ages |${\gt}8 \, \rm {Gyr}$|⁠), i.e. assuming the only cluster mass loss was due to stellar evolution. On average (individual galaxies may of course differ), for galaxies with |$M_\ast \lt 10^{10} \, \rm {M}_{\odot }$| the z = 0 GC metallicity distributions (solid black lines) are a reasonable reflection of the initial distributions, but typically a factor 3–4 lower due to dynamical mass loss. However, at higher masses there is a clear dependence of GC ‘survival’ on metallicity (at least above the adopted GC luminosity limit), with surviving high-metallicity (⁠|$[\rm {Fe}/\rm {H} ]\sim 0$|⁠) clusters representing a smaller fraction (about one-tenth) of the initial population than at lower metallicities. This effect also appears to increase in higher-mass galaxies. We will return to this point in Section 3.2 when discussing the specific frequency–metallicity relation.

In Fig. 2 the very massive galaxies (⁠|${\gt}10^{11} \, \rm {M}_{\odot }$|⁠) appear to have too many low metallicity clusters. In Fig. 3 we show this is due to the small ACS Fornax Cluster Survey field of view (202 × 202 arcsec), combined with the radial gradients in GC metallicity. Adopting a similar radius limit (projected radius of |$20 \, \rm {kpc}$|⁠), we find excellent agreement with the metallicity distributions of Fornax galaxies. Potentially, at |$[\rm {Fe}/\rm {H} ]\approx -3$| the number of simulated GCs are still over-predicted by a factor of ≈3 (which is similarly suggested in the 10.5 < log₁₀(M_*/M_⊙) < 11 panel of Fig. 2). Such particles form in very low mass galaxies (i.e. |$M_\ast \lesssim 10^6 \, \rm {M}_{\odot }$|⁠, as inferred in fig. 9 of Kruijssen et al. 2019a) which are poorly resolved in the simulation (with initial baryonic particle masses of |$2.25\times 10^5 \, \rm {M}_{\odot }$|⁠). Thus an over-abundance of extremely low-metallicity GCs could stem from the GC formation model performing incorrectly in poorly-resolved galaxies (e.g. overestimating |$M_{\rm c,\ast }$|⁠).

3.2 Specific frequency–metallicity relationship

In Fig. 4 we compare the relationship between specific frequency (i.e. number of GCs per unit galaxy mass; T_N = N_GC/(M_*/M_⊙)) and metallicity in the E-MOSAICS galaxies. The simulations predict a strong dependence of T_N on metallicity, with more GCs per unit mass at lower metallicities, as found for observed galaxies (Harris & Harris 2002; Lamers et al. 2017). Across all stellar mass ranges, the slope of the T_N-|$[\rm {Fe}/\rm {H} ]$| relations for simulated galaxies (at least for GC ages |${\gt}8 \, \rm {Gyr}$|⁠) is in reasonable agreement with that found for NGC 5128 and M31 (Lamers et al. 2017). The simulations predict slightly flatter slopes for the T_N-|$[\rm {Fe}/\rm {H} ]$| relation at higher galaxy masses, though we caution that this result may be strongly affected by the under-disruption of GCs in the simulations. As the EAGLE model does not resolve the cold, dense phase of the ISM, mass loss through tidal shocks by dense substructure is under-estimated in the E-MOSAICS simulations, particularly at higher metallicities (see Pfeffer et al. 2018; Kruijssen et al. 2019a, for further details and discussion).

In galaxies with stellar masses |$M_\ast \lt 10^{10} \, \rm {M}_{\odot }$|⁠, the T_N-|$[\rm {Fe}/\rm {H} ]$| relation (solid line) generally follows the initial relation (dotted line); i.e. the T_N-|$[\rm {Fe}/\rm {H} ]$| relation is largely set by cluster formation alone. In the E-MOSAICS model, this increase of T_N towards low metallicities at the time of GC formation is driven by the CFE (fraction of star formation in bound clusters), which tends to increase at lower metallicities (Pfeffer et al. 2018). This increase is a result of the metallicity-dependent density threshold for star formation implemented in EAGLE, which accounts for the transition from a warm, photoionized interstellar gas to a cold, molecular phase, that is expected to occur at lower densities/pressures at higher metallicities (Schaye 2004).

Towards higher galaxy masses, the initial relations show an increase in T_N at |$[\rm {Fe}/\rm {H} ]\sim 0$|⁠. Cluster mass loss then begins to play a more significant role in setting GC numbers, particularly at |$[\rm {Fe}/\rm {H} ]\gtrsim -1$| (as noted in Section 3.1; see Kruijssen 2015), i.e. cluster mass loss preferentially affects higher-metallicity GCs in higher-mass galaxies. Both processes are due to an increase in higher density/pressure gas at higher galaxy masses, resulting in both increased cluster formation (higher CFE) and cluster disruption (i.e. the ‘cruel cradle effect’; Kruijssen et al. 2012). We note that in other GC formation models the specific frequency–metallicity relation may have different origins. For example, in the Choksi et al. (2018) model the specific frequency–metallicity relation is set by cluster formation (GC disruption in the model is environmentally independent), while they find a good match to the observed slope of the relation.

Observations of the T_N-|$[\rm {Fe}/\rm {H} ]$| relation are generally restricted to galaxy haloes and exclude the galaxy discs (e.g. |$R \gt 6 \, \rm {kpc}$| for NGC 5128 and |$R \gt 12 \, \rm {kpc}$| for M31; Harris & Harris 2002; Lamers et al. 2017). Therefore in Fig. 5 we compare the relations as a function of (3D) radius within the simulated galaxies. To account for different sizes of galaxies, we compare four radial ranges scaled by the 3D stellar half-mass radius (0–1, 1–2, 2–5, and 5–10r_h,*).

For galaxy stellar masses |$M_\ast \lt 10^9 \, \rm {M}_{\odot }$| we do not find strong differences in the T_N-|$[\rm {Fe}/\rm {H} ]$| relation as a function of radius. In the range 10⁹-|$10^{10.5} \, \rm {M}_{\odot }$| there are clear radial gradients in T_N at |$[\rm {Fe}/\rm {H} ]\gtrsim -1$|⁠, such that T_N is higher at smaller radii. This is driven by radial gas pressure gradients in galaxies, which results in higher CFE at smaller radii (see fig. 4 in Pfeffer et al. 2018), and an ex situ origin for GCs at larger radii through the accretion of lower mass galaxies (fig. 8 in Reina-Campos et al. 2022a). At radii >5r_h,* (i.e. the ‘halo’ of the galaxies) the T_N-|$[\rm {Fe}/\rm {H} ]$| relations are in good agreement with the observations from Lamers et al. (2017) due to this (largely) ex situ origin, despite the under-disruption of GCs at smaller radii.

For galaxy stellar masses |$M_\ast \gt 10^{10.5} \, \rm {M}_{\odot }$| we find that T_N at high metallicities (⁠|$[\rm {Fe}/\rm {H} ]\gtrsim -1$|⁠) is increased at all radii. This is due to the increased importance of mergers, particularly major mergers, in the assembly of massive galaxies (Qu et al. 2017). As a result, an ex situ origin for GCs becomes dominant in massive galaxies at nearly all radii (Reina-Campos et al. 2022a; Trujillo-Gomez et al. 2022) and metallicities (Section 4.1). The merging of massive galaxies (with elevated T_N at high |$[\rm {Fe}/\rm {H} ]$| and small radii) then results in results in elevated T_N at high |$[\rm {Fe}/\rm {H} ]$| and large radii in the merger descendant.

Interestingly, comparison of the T_N-|$[\rm {Fe}/\rm {H} ]$| relations (Figs 4 and 5) and GC metallicity distributions (Figs 2 and 3) for the highest mass galaxies (⁠|$M_\ast \gt 10^{11} \, \rm {M}_{\odot }$|⁠) suggests that the T_N-|$[\rm {Fe}/\rm {H} ]$| relation may differ in Fornax cluster galaxies. While the simulated T_N-|$[\rm {Fe}/\rm {H} ]$| relation is flatter than observed in M31 and NGC 5128 (and similarly, the number of high-metallicity GCs is larger than observed in M31; Fig. 2), we find good agreement with the metallicity distributions for Fornax cluster galaxies (Fig. 3). The origin of this difference in the observations is not clear, but plausibly could be due to differences in the formation histories of cluster and field galaxies (such as differing galaxy merger or star-formation quenching histories). Ideally, future observations should test for an environmental dependence in the T_N-|$[\rm {Fe}/\rm {H} ]$| relations of galaxies (which might be possible in the Virgo cluster using observations of individual stars from future extremely large telescopes, e.g. Deep et al. 2011).

3.3 Effect of cluster formation model

A benefit of the E-MOSAICS approach to modelling galaxy and GC formation is that multiple star cluster formation models can be tested simultaneously in parallel within the same simulation. In this section we test the effect of the different assumptions about the formation physics of GCs on the resulting GC metallicity distributions. We consider four formation models (i.e. two settings for each of the CFE and |$M_{\rm c,\ast }$| models): varying CFE and |$M_{\rm c,\ast }$| (fiducial model), constant CFE and power-law initial mass function (‘constant formation’ model), varying CFE and power-law initial mass function (‘CFE-only’ model), and constant CFE and varying |$M_{\rm c,\ast }$| (‘|$M_{\rm c,\ast }$|-only’ model).

In Fig. 6, we compare the metallicity distributions for the four formation models as a function of galaxy mass (as in Fig. 2). As discussed by Bastian et al. (2020), the number of GCs predicted varies significantly between the formation models, with CFE-only and constant formation models (i.e. variations without an upper mass truncation) predicting the most GCs and the |$M_{\rm c,\ast }$|-only model predicting the fewest. The shapes of the GC metallicity distributions also vary between the different models. At galaxy masses |$M_\ast \lt 10^9 \, \rm {M}_{\odot }$|⁠, the three alternative formation models peak at higher metallicities than the fiducial model, while at |$M_\ast \gtrsim 10^{10} \, \rm {M}_{\odot }$| the constant-formation model peaks at lower metallicities than the fiducial model.

We quantify this further in Fig. 7, by comparing the median GC metallicity as a function of galaxy mass with observed galaxies. We find the fiducial model (top left-hand panel) agrees well with the median GC metallicity in observed galaxies, other than at |$M_\ast \sim 10^{10.5} \, \rm {M}_{\odot }$| where the simulated galaxies produce too many high-metallicity GCs. At very large masses (⁠|$M_\ast \gt 10^{11} \, \rm {M}_{\odot }$|⁠) the simulated galaxies have median GC metallicities much lower than the Fornax cluster galaxies. As discussed in Section 3.1, this is related to limited ACS Fornax Cluster Survey field of view. In Appendix  A we show the predicted median metallicities when adopting a |$20 \, \rm {kpc}$| radius limit. With this limit the median GC metallicities of massive galaxies are in good agreement with Fornax/Virgo clusters and SLUGGS survey galaxies, but over-predict the metallicities relative to field galaxies (Milky Way, M31, M81).

The comparison of the four different formation models in Fig. 7 reveals how the CFE and |$M_{\rm c,\ast }$| models affect the GC metallicity–galaxy mass relation. The models with an environmentally-dependent CFE (fiducial, CFE-only) produce relations with slopes in reasonable agreement with the observed galaxies, while models with a constant CFE (constant, |$M_{\rm c,\ast }$|-only) produce relations with slopes that are too flat. This is because star formation generally occurs at higher pressures in higher mass galaxies (at least at the peak of the star-formation rate), thus resulting in higher CFE in the environmentally-dependent model (though this will vary significantly depending on galaxy formation history, for further discussion see Pfeffer et al. 2018). Models without an upper GC mass function truncation (constant, CFE-only) generally have too high GC metallicities for low-mass galaxies (⁠|$M_\ast \lt 10^{10} \, \rm {M}_{\odot }$|⁠), while models with an environmentally-dependent |$M_{\rm c,\ast }$| are in better agreement with observed galaxies at the same masses (i.e. the |$M_{\rm c,\ast }$| model may suppress GC formation at high metallicities). Thus, the CFE model influences the slope of the GC metallicity–galaxy mass relation, while the upper mass function truncation (or lack of) influences the normalization.

Comparison of the GC formation models in Fig. 6 also suggests that the underabundance of low metallicity GCs (⁠|$[\rm {Fe}/\rm {H} ]\lt -1$|⁠) in galaxies with masses 9.5 < log (M_*/M_⊙) < 10 in the fiducial model (noted in Section 3.1) could be due to a suppression of GC numbers by |$M_{\rm c,\ast }$|⁠. Comparing the fiducial and CFE-only models (i.e. with and without an |$M_{\rm c,\ast }$| model), the CFE-only model appears in better agreement with the Fornax cluster galaxies. However, we note that overall the upper truncation masses of GC systems in the fiducial model are in very good agreement with Virgo cluster galaxies, including this mass range (Hughes et al. 2022). Additionally, the CFE-only model over-predicts low metallicity GC numbers in higher mass galaxies (⁠|$M_\ast \gt 10^{10.5} \, \rm {M}_{\odot }$|⁠). In the future, detailed comparisons of the upper truncation mass of GC systems as a function of metallicity and galaxy mass between observed and simulated galaxies may help to resolve the issue. The problem could also stem from the EAGLE model (i.e. too few low metallicity star particles, within which the GCs form) as galaxies in the ‘Recalibrated’ EAGLE model have median stellar metallicities ∼0.2–0.3 dex higher than observed within this mass range (which itself could also be due to uncertainty in nucleosynthetic yields or absence of metal mixing between SPH particles Schaye et al. 2015). Alternatively, our GC luminosity selection (Section 2.2) may be too bright in comparison to the observed galaxies in this mass range.

3.4 GC metallicity distribution bimodality

We now test whether the E-MOSAICS GC formation model produces bimodal GC metallicity distributions in a way that is consistent with observations. We use the Gaussian mixture modelling (GMM) algorithm from Muratov & Gnedin (2010) to test whether a unimodal or bimodal distribution is preferred for the GC metallicity distribution of each galaxy with at least 30 GCs. Following Muratov & Gnedin (2010), we consider a distribution to be bimodal if the bimodal solution is preferred with a probability p > 0.9, the distribution has a negative kurtosis and the two peaks are separated by a relative distance D > 2 (defined as D = ‖μ₁ − μ₂‖/[(σ₁² + σ₂²)/2]^1/2, where μ and σ are the mean and standard deviation of a Gaussian distribution, respectively).

In Fig. 8 we compare the outcomes of the bimodality tests, showing the metallicities of the inferred ‘blue’ (low metallicity, top left-hand panel) and ‘red’ (high metallicity, top right-hand panel) peaks and ‘red’ GC fraction for bimodal galaxies (bottom left-hand panel), as well as the fraction of bimodal galaxies as a function of galaxy mass (bottom right-hand panel). We show the results for both GC age limits of |${\gt}8 \, \rm {Gyr}$| (standard sample, black lines) and |${\gt}2 \, \rm {Gyr}$| (including younger disc GCs, orange lines). Comparing the GC age |${\gt}8 \, \rm {Gyr}$| samples (the most relevant comparison for cluster galaxies given they are generally old due to early star formation quenching, e.g. Gallazzi et al. 2021), we find very good agreement with the Virgo cluster results from Peng et al. (2006) for the trend of mean blue and red GC peak metallicities with galaxy mass (top panels), showing a trend of higher-metallicity peaks in higher-mass galaxies. The drop or flattening in peak metallicities at high galaxy masses, relative to the Virgo cluster trend, is a result of the absence of a radius limit in the GC selection (see Section 3.3 and Appendix  A). A similar flattening in the location of the peak metallicities was found for galaxies with |$M_\ast \gtrsim 10^{11} \, \rm {M}_{\odot }$| (halo masses |${\gtrsim}10^{13} \, \rm {M}_{\odot }$|⁠) by Choksi et al. (2018). With an age limit of |${\gt}2 \, \rm {Gyr}$| the peak metallicities are ≈0.1–0.2 dex higher, with the largest difference for the metal-rich peak in low-mass galaxies. This is a result of galaxy downsizing, where lower-mass galaxies form more of their mass at later times (Bower et al. 1992; Gallazzi et al. 2005).

The bottom left-hand panel of Fig. 8 shows the fraction of GCs in the red (high metallicity) peak for each bimodal galaxy. For galaxy stellar masses |${\gt} 10^{9.5} \, \rm {M}_{\odot }$| and a |${\gt}8 \, \rm {Gyr}$| age limit we find good agreement with Virgo cluster galaxies (Peng et al. 2006), with a red fraction that increases with galaxy mass and plateauing at ≈0.55. With a |${\gt}2 \, \rm {Gyr}$| age limit the red fraction is only slightly higher than for a |${\gt}8 \, \rm {Gyr}$| limit. At lower galaxy masses the red fraction appears to increase (for ages |${\gt}8 \, \rm {Gyr}$|⁠), unlike the observed fraction. For a number of these low-mass galaxies the ‘red’ peak occurs at |$[\rm {Fe}/\rm {H} ]\sim -1.4$|⁠, which is similar to the blue peak in most other galaxies. The occurrence of a second ‘blue’ peak at lower metallicities in these galaxies could just be due to randomness where the number of GCs is low (i.e. ∼30 GCs in total). Indeed, for galaxies with a blue mean |$[\rm {Fe}/\rm {H} ]\sim -2.5$| the typical number of GCs in the blue peak is only ∼10 (with a range of 3–20). However, we also caution that at such masses (⁠|${\sim}10^9 \, \rm {M}_{\odot }$|⁠) we are comparing a biased sample of galaxies due to our requirement for a minimum number of GCs (30) to perform the bimodality tests. The fraction of galaxies with at least 30 GCs is shown as a thin dashed line in the figure and is <10 per cent at |$M_\ast = 10^9 \, \rm {M}_{\odot }$|⁠. At similar masses, the median GC metallicity is |$[\rm {Fe}/\rm {H} ]\approx -1.5$| (Fig. 7), thus the selection on GC numbers leads to a highly-biased sample of low-mass galaxies with a significant number of metal-rich GCs. For completeness, the dotted lines in the figure show the results of the GMM tests for all galaxies, rather than only bimodal galaxies, though considering blue and red ‘peaks’ for a unimodal distribution is not meaningful and simply reflects that the distribution is not a perfect Gaussian. Comparing the GMM tests for all galaxies does not significantly change the peak metallicities, but does increase the red fractions to >50 per cent.

As an alternative, more general, comparison, in Fig. 9 we compare the fraction of red (high metallicity) GCs using a fixed split between the populations at |$[\rm {Fe}/\rm {H} ]= -1$| (a similar approach was used by Choksi & Gnedin 2019b for comparing non-bimodal galaxies). Though using fixed metallicity split ignores the increasing mean GC metallicity with galaxy mass, it can be applied to all galaxies, rather than only those with bimodal GC populations. Similarly, it avoids issues where the ‘red’ GC peak for some galaxies occurs at similar metallicities to typical blue GC peaks. For a GC age limit |${\gt}8 \, \rm {Gyr}$| (grey points, solid black line), the simulations show good agreement with observed galaxies at galaxy stellar masses |${\lesssim}10^{9.5} \, \rm {M}_{\odot }$| and |${\gtrsim}10^{10.5} \, \rm {M}_{\odot }$|⁠. At |$M_\ast \approx 10^{10} \, \rm {M}_{\odot }$|⁠, simulated galaxies tend to have elevated red GC fractions by ≈0.2, consistent with the elevated median GC metallicities at similar masses (Fig. 7). With a GC age limit of |${\gt}2 \, \rm {Gyr}$| (solid orange line), we find elevated red GC fractions compared to using a |${\gt}8 \, \rm {Gyr}$| limit, with a typical difference in the fractions ranging from 0.1 to 0.4 depending on galaxy mass (peaking at |${\approx}10^{9.5} \, \rm {M}_{\odot }$|⁠). A difference of ≈0.34 is also found for the LMC, with a smaller difference (≈0.07) for the SMC (red squares connected by vertical lines). The differences between the simulation trends in Fig. 9 and the bottom right-hand panel of Fig. 8 (e.g. the increasing red fraction for low-mass galaxies in Fig. 8) can be explained by galaxies with both low and high red fractions tending not to have bimodal GC populations.

In the bottom right-hand panel of Fig. 8 we show the fraction of galaxies with bimodal GC metallicity distributions as a function of galaxy mass. We find that 37 ± 2 per cent of simulated GC distributions are bimodal (for an age limit |${\gt} 8 \, \rm {Gyr}$|⁠) and the fraction is slightly higher (45 ± 7 per cent) for higher mass galaxies (⁠|$M_\ast \gt 10^{10.5} \, \rm {M}_{\odot }$|⁠). For an age limit |${\gt} 2 \, \rm {Gyr}$|⁠, the bimodal fraction is significantly lower (≈10 per cent) at masses |${\sim}10^{10} \, \rm {M}_{\odot }$|⁠. This difference is largely caused by galaxies failing the kurtosis test (i.e. having a positive kurtosis) due to the increase in metal-rich GC numbers. For comparison we also show the combined bimodal fractions from the Fornax cluster (Fahrion et al. 2020b), SLUGGS survey (Usher et al. 2012), and Local Group galaxies (MW, M31, M81). The total bimodal fraction for the observed galaxies is |$44^{+10}_{-9}$| per cent (12/27 galaxies). Overall, we find the predicted bimodal fractions from the simulations are in very good agreement with observed fractions.

To understand the origin of bimodal GC distributions in the simulated galaxies, we compare the metallicity distributions of two galaxies with clear bimodal distributions in Fig. 10. The galaxies were selected to be central galaxies in Milky Way-mass haloes (⁠|$M_{200} \approx 10^{12} \, \rm {M}_{\odot }$|⁠). Despite relatively similar GC metallicity distributions for both galaxies, comparing the z = 0 (solid line) and initial (dash–dotted line) distributions in each case shows a very different origin for bimodality.

In Group 86 (bottom panel) the bimodal GC distribution appears to be imprinted at cluster formation, with the initial distribution appearing as a scaled-up version of the z = 0 distribution. However, this is not reflected in the stellar metallicity distribution (for the same |${\gt}8 \, \rm {Gyr}$| age limit), which is unimodal and peaks at a similar metallicity to the ‘red’ GC peak (⁠|$[\rm {Fe}/\rm {H} ]\approx -0.5$|⁠). Interestingly, when also considering younger disc GCs (ages |${\gt} 2 \, \rm {Gyr}$|⁠) the GC metallicity distribution of this galaxy appears to be trimodal (with a third peak at |$[\rm {Fe}/\rm {H} ]\approx 0.2$|⁠). For this galaxy, the majority of the ‘blue’ GC peak is in fact made up of in situ clusters (unlike for Group 85), though there is a clear gradient in ex situ fraction as a function of metallicity, which we will discuss further in Section 4.1.

In Group 85 (top panel), neither the initial GC or stellar metallicity distributions are bimodal. Instead, the GC distribution evolves to be bimodal through dynamical mass loss. We quantify this further in Fig. 11, by comparing the ratio of the final (z = 0) to initial number of GCs as a function of metallicity for central galaxies within |$M_{200} \approx 10^{12} \, \rm {M}_{\odot }$| haloes. For Group 86, GC mass loss is similar to the typical mass loss for such galaxies, being slightly more effective at high metallicities (⁠|$[\rm {Fe}/\rm {H} ]\sim 0$|⁠) than at lower metallicities (⁠|$[\rm {Fe}/\rm {H} ]\sim -2$|⁠; a result of the more disruptive, higher-density environment in which metal-rich GCs are formed; see Pfeffer et al. 2018 and Kruijssen et al. 2019a for further discussion within the context of the E-MOSAICS simulations). However, for Group 85 cluster mass loss is particularly effective at |$[\rm {Fe}/\rm {H} ]\sim -1$| (i.e. at the trough in the z = 0 distribution), which results in the evolution into a bimodal metallicity distribution.

We investigate the origin of the bimodal initial GC distribution for Group 86, and lack of one in the initial distribution for Group 85, in Fig. 12 by comparing the distributions of the cluster formation properties CFE (left-hand panel) and |$M_{\rm c,\ast }$| (right-hand panel) as a function of metallicity. Both galaxies show a decline in CFE from |$[\rm {Fe}/\rm {H} ]= -3$| to −1, meaning GCs are least likely to form at |$[\rm {Fe}/\rm {H} ]\approx -1$|⁠. This decline is driven by the decreasing density threshold for star formation (see fig. 3 in Pfeffer et al. 2018) which models the transition from a warm to cold ISM that is expected to occur at lower densities for higher metallicities (Schaye 2004), meaning at very low metallicities the fraction of star formation in clusters (CFE) is high. The CFE in both galaxies then experiences an increase at higher metallicities as a result of high natal gas pressures at small galactic radii (which also tends to result in GCs with elevated α-abundances; Hughes et al. 2020). However, the initially bimodal GC distribution in Group 86 is not completely caused by a change in the CFE (at least in this case), as the typical CFE (≈10 per cent) is relatively similar across the range |$-1.5 \lt [\rm {Fe}/\rm {H} ]\lt -0.5$| with only a small decrease at |$[\rm {Fe}/\rm {H} ]\approx -1.1$|⁠. Instead, it appears to be further driven by an evolution in |$M_{\rm c,\ast }$| (right-hand panel in the figure, which itself also depends on the CFE), such that fewer GCs are formed at high enough masses that they pass the luminosity selection at z = 0. We find that the typical |$M_{\rm c,\ast }$| for Group 86 falls below |$10^5 \, \rm {M}_{\odot }$| at |$[\rm {Fe}/\rm {H} ]\approx -1.1$|⁠, i.e. at the same metallicity where the trough in the initial (and indeed final) GC metallicity distribution is found for this galaxy. This is not the case for Group 85, where the typical |$M_{\rm c,\ast }$| remains relatively similar (⁠|${\approx}10^{5.5} \, \rm {M}_{\odot }$|⁠) at metallicities |$[\rm {Fe}/\rm {H} ]\lt -0.5$|⁠. Therefore, for some galaxies, variations in cluster formation properties (CFE and |$M_{\rm c,\ast }$|⁠) as a function of metallicity can result in initially-bimodal GC metallicity distributions.

4 DISCUSSION

4.1 Origin of bimodal GC metallicity distributions

Most previous works have focussed on GC formation as the origin of bimodal metallicity distributions or otherwise do not separately consider the effects of cluster formation and disruption. In general, bimodality has been suggested to occur through discrete episodes of GC formation, e.g. via galaxy mergers (e.g. Zepf & Ashman 1993; Muratov & Gnedin 2010) or truncating the formation of metal-poor GCs below some redshift (e.g. Forbes et al. 1997; Beasley et al. 2002; Strader et al. 2005; Griffen et al. 2010). Models where GC disruption is environmentally independent have nevertheless been successful in reproducing the observed peak metallicities of red and blue GC populations (Choksi et al. 2018).

In the E-MOSAICS model GCs may form in a galaxy where star formation is sufficiently intense, and thus there is no explicit truncation in GC formation times (other than that which occurs ‘naturally’ during the formation of galaxies, such as star formation quenched by stellar/AGN feedback or gas exhaustion, which will therefore vary significantly from galaxy to galaxy depending on their formation histories). As shown in Section 3.4, in the E-MOSAICS simulations GC metallicity bimodality arises in different ways: through variations in cluster formation (efficiency and upper mass truncation) as a function of metallicity creating an initially bimodal distribution, or cluster disruption creating a bimodal distribution from an initially unimodal distribution.

To determine whether cluster formation or disruption is more important in the origin of bimodal metallicity distributions in the E-MOSAICS model, in Fig. 13 we show the results of performing GMM tests on the initial GC metallicity distributions for galaxies in Fig. 8. Unlike for the z = 0 bimodal fractions which are nearly constant at all galaxy masses (Fig. 8), the initial bimodal fractions strongly depend on galaxy mass. For |$M_\ast \gtrsim 10^{9.5} \, \rm {M}_{\odot }$| (where >50 per cent of the galaxy population is sampled) the initial bimodal fractions decline with increasing galaxy mass (≈60 to 0 per cent from 10^9.5 to |$10^{10.5} \, \rm {M}_{\odot }$|⁠). In fact for galaxies with stellar masses |$M_\ast \gt 10^{10.5} \, \rm {M}_{\odot }$|⁠, no galaxies have initially bimodal GC metallicity distributions, despite ≈40 per cent of such systems being bimodal at z = 0.

For |$M_\ast \gtrsim 10^{10} \, \rm {M}_{\odot }$| (where the initial bimodal fraction falls below the z = 0 fraction), Fig. 13 shows that cluster mass loss is therefore the most important factor for determining whether a GC metallicity distribution will become bimodal in the E-MOSAICS model. A similar point was made by Kruijssen (2015) who argued that GC destruction, which sets the specific frequency–metallicity relation (Section 3.2), is responsible for observed GC metallicity distributions. However, for galaxies with masses |$M_\ast \lt 10^{10} \, \rm {M}_{\odot }$| the inverse appears to be the case: More galaxies (47 per cent) have initially bimodal GC distributions than at z = 0 (37 per cent). Thus in low-mass galaxies GC evolution decreases the fraction of bimodal galaxies.

In principle, determining whether formation or disruption resulted in the GC metallicity distribution of observed galaxies could be tested by surveying the field stars from destroyed GCs. In fact, the multiple stellar populations observed in Galactic GCs (light-element abundance spreads that appear to be unique to GCs; see reviews by Gratton et al. 2012; Charbonnel 2016; Bastian & Lardo 2018) potentially provide such a method. By searching for their unique chemical fingerprints, former GC stars have been identified in the stellar halo (Martell & Grebel 2010; Martell et al. 2016; Koch, Grebel & Martell 2019; Horta et al. 2021b) and bulge (Schiavon et al. 2017) of the Milky Way.

Interestingly, in the inner Galaxy the nitrogen-rich (potentially former GC) stars peak at a metallicity |$[\rm {Fe}/\rm {H} ]\approx -1$| (Schiavon et al. 2017), unlike the Galactic GC system which has a trough at a similar metallicity (see fig. 9 in Schiavon et al. 2017). Thus the N-rich stars do not appear to originate from surviving GCs, but instead could be the remnants of destroyed GCs. This suggests that the bimodal metallicity distribution of Milky Way’s GC system may be a result of GC disruption, rather than formation, in line with the expectations for massive galaxies in the E-MOSAICS model.

In the context of hierarchical galaxy formation, it has been suggested that metal-rich and metal-poor GCs represent in situ and ex situ formation (e.g. Côté et al. 1998; Hilker et al. 1999; Forbes & Remus 2018). We test this idea in Fig. 14 by comparing the average ex situ GC fractions as a function of metallicity and galaxy mass. Star particles (and their associated GCs) are defined as in situ or ex situ based on the subhalo the parent gas particle was bound to in the snapshot immediately prior to star formation (i.e. the time when the gas particle was converted into a star particle). Particles formed in subhaloes within the main progenitor branch are defined as in situ formation.

The metallicity below which ex situ GCs become dominant (fraction >50 per cent) shows a strong correlation with galaxy mass, from |$[\rm {Fe}/\rm {H} ]\approx -2.3$| in low-mass galaxies (⁠|$M_\ast \lt 10^{8.5} \, \rm {M}_{\odot }$|⁠) to |$[\rm {Fe}/\rm {H} ]\approx 0$| in massive galaxies (⁠|$M_\ast \gt 10^{11} \, \rm {M}_{\odot }$|⁠). We show this directly in Fig. 15. This is a consequence of the increasing contribution of mergers to the growth of galaxies at higher masses (e.g. Qu et al. 2017; Clauwens et al. 2018; Pillepich et al. 2018; Tacchella et al. 2019; Davison et al. 2020). As shown by Choksi & Gnedin (2019b), the fraction of ex situ GCs reasonably follows the increasing ex situ field star fraction with increasing galaxy mass. Thus in massive galaxies (⁠|$M_\ast \gt 10^{11} \, \rm {M}_{\odot }$|⁠) GC bimodality is not directly related to in/ex situ GC formation as the majority of GCs are formed ex situ, even in the red peak. Rather, the GC metallicities are related to the time of formation and the galaxy mass-metallicity relation (Muratov & Gnedin 2010; Kruijssen 2015, 2019). For galaxies with masses |${\lesssim}5 \times 10^9 \, \rm {M}_{\odot }$|⁠, ex situ GCs only become dominant at metallicities below the typical blue peak metallicity. Only in galaxies with stellar masses of |${\approx}10^{10.5} \, \rm {M}_{\odot }$| will the blue and red peaks (generally) correspond to ex situ and in situ GC formation. Naturally, this will also vary from galaxy to galaxy, as shown in Fig. 10 where ex situ GCs dominate the blue peak for Group 85, but in situ GCs dominate the blue peak for Group 86. A correlation between the accreted fractions of metal-poor/rich GCs and galaxy mass was similarly found by Choksi & Gnedin (2019b), despite the very different GC models, and thus appears to be a general feature of cosmological GC formation models based on ‘young star cluster’ formation (see next section).

4.2 GC metallicity distributions and the origin of GCs

The metallicity distributions of GC systems, bimodal distributions in particular, have been seen as critical in understanding the origin of GCs, with a number of different scenarios proposed to explain their origin. Some works have investigated a direct connection between GC formation and gas-rich galaxy mergers (Ashman & Zepf 1992; Zepf & Ashman 1993; Beasley et al. 2002; Muratov & Gnedin 2010), inspired by the the young star clusters with GC-like properties observed in merging and post-merger galaxies (e.g. Holtzman et al. 1992; Whitmore et al. 1993; Whitmore & Schweizer 1995). In this scenario, many mergers of low-mass galaxies and few mergers of high-mass galaxies thus create the metal-poor and metal-rich GC distribution peaks, respectively. Such models can reproduce many scaling relations of GC populations (Li & Gnedin 2014; Choksi et al. 2018). However, young massive star clusters have been observed forming in a variety of environments (e.g. from regular spiral galaxies to star-burst dwarf galaxies), not just merging and interacting galaxies (O’Connell, Gallagher & Hunter 1994; Larsen & Richtler 1999; Hunter et al. 2000; Billett, Hunter & Elmegreen 2002; Mora et al. 2009; Chandar et al. 2010; Annibali et al. 2011). Though galaxy mergers can help to drive conditions favourable for star cluster formation (higher gas pressures and star formation rates, e.g. Bekki et al. 2002; Bournaud, Duc & Emsellem 2008; Lahén et al. 2020; Li et al. 2022), massive star cluster formation may still occur in the absence of galaxy mergers. Mergers may also promote the survival of star clusters by ejecting them from the disruptive environment of star-forming discs (Kravtsov & Gnedin 2005; Kruijssen 2015), and in this way merger-based models may have GC survival included in the formation model. Updated versions of the Muratov & Gnedin (2010) semi-analytic model (which model galaxies and GC populations by applying analytic scaling relations to dark matter-only cosmological simulations) instead consider GC formation during high halo accretion rates, rather than just during mergers (Choksi et al. 2018; Choksi & Gnedin 2019b).

Other works have suggested separate formation mechanisms for metal-poor and metal-rich GCs. While metal-rich GCs may form in mergers or otherwise with the bulk of the field stars, metal-poor GCs have been suggested to form in low-mass dark matter haloes in the early Universe (Peebles 1984; Rosenblatt et al. 1988; Ricotti, Parry & Gnedin 2016), whose formation is possibly truncated (e.g. Strader et al. 2005; Moore et al. 2006; Griffen et al. 2010) or caused (Cen 2001) by reionization. As already discussed in previous works, the scenario has a number of fundamental issues. Isolated GCs, residing in dark matter haloes yet to be accreted into a larger system, have not been observed (di Tullio Zinn & Zinn 2015; Mackey, Beasley & Leaman 2016). Some GCs (i.e. those which always resided in weak tidal fields) might be expected to retain their dark matter halo, however GCs at large distances from their host galaxy are consistent having no dark matter halo (Baumgardt et al. 2009; Lane et al. 2010; Conroy, Loeb & Spergel 2011). Critically, to explain GC metallicities such GCs must self-enrich through supernovae feedback, which is inconsistent with the very small iron spreads (<0.1 dex; Carretta et al. 2009) in most² GCs (Peebles 1984; Rosenblatt et al. 1988).

In this work, we show that a GC formation and evolution model which reproduces the observed scaling relations of young star cluster populations (using the fiducial E-MOSAICS model; Pfeffer et al. 2019) also reproduces many properties of GC metallicity distributions when applied to galaxy formation in a cosmological context. The model does not rely on separate formation mechanisms or truncated formation epochs for metal-poor and metal-rich GCs, nor preferential formation in galaxy mergers (which generally are not a major contribution to GC formation in the E-MOSAICS model; Keller et al. 2020). Instead, variations in GC metallicity distributions between galaxies may arise due to local environmental variations in GC formation properties and their subsequent evolution over time.

A ‘young star cluster’ scenario for the origin of GCs can also explain many other properties of GC populations in galaxies, including the masses of GCs (Kravtsov & Gnedin 2005; Kruijssen 2015; Li et al. 2017; Ma et al. 2020), the fraction of stars contained in GCs (Choksi et al. 2018; Bastian et al. 2020), the typical old ages of GCs (Reina-Campos et al. 2019), the age–metallicity relations of GC systems (Kruijssen et al. 2019a; Li & Gnedin 2019; Horta et al. 2021a), the radial distributions of GC systems (Reina-Campos et al. 2022a), the high-mass truncation of GC mass functions (Hughes et al. 2022), the redder colours of more luminous GCs (‘blue tilt’; Usher et al. 2018; Choksi & Gnedin 2019a; Kruijssen 2019), and the UV luminosity function of proto-GCs at high redshifts (Bouwens et al. 2021). The uncertainties that remain in this scenario are generally related to the difficulties in modelling both the small scales of GCs (∼parsecs) and the large scales of galaxy populations (∼megaparsecs) or particular numerical implementations (e.g. such as the overabundance of low-mass GCs in the E-MOSAICS simulations; Pfeffer et al. 2018). Considered together, this body of work suggests a common origin for GCs and young star clusters.

5 SUMMARY

In this work, we investigate the GC metallicity distributions predicted by the E-MOSAICS simulations. The main aim of the work is to test whether a common formation scenario for young star clusters and GCs can explain the observed GC metallicity distributions. We find that:

• The predicted GC metallicity distributions from the fiducial E-MOSAICS model generally agree well with the observed distributions as a function of galaxy mass (Section 3.1). In particular, the distributions for high-mass galaxies (⁠|$M_\ast \gt 10^{11} \, \rm {M}_{\odot }$|⁠) are in very good agreement with observed distributions once the field of view is accounted for (Fig. 3). However, simulated Milky Way-mass galaxies (⁠|$10.5 \lt \log (M_\ast / \rm {M}_{\odot }) \lt 11$|⁠) often have too many metal-rich GCs compared to observed galaxies (as already discussed in detail in Kruijssen et al. 2019a), which is suggested to be due to inadequate disruption of GCs through lack of substructure in galactic discs (see Section 2.1 for a brief summary). The number of extremely low metallicity GCs (⁠|$[\rm {Fe}/\rm {H} ]\lesssim -2.5$|⁠) are also overpredicted (by a factor ≈3), which may be a result of the resolution limits of the simulation.

• The predicted relationship between GC specific frequency (number of GCs per unit mass) and metallicity shows a slope that is consistent with observed trends (Section 3.2). In the simulations, this trend emerges due to an increased CFE at low metallicities and preferential disruption of high metallicity GCs in higher mass galaxies. However, the under-disruption of high-metallicity GCs at small galactocentric radii is also evident in the specific frequency–metallicity relations.

• Comparing the variations of the E-MOSAICS GC formation model (Section 3.3), we find that models with a constant CFE fail to reproduce the correlation between median GC system metallicity and galaxy mass, instead predicting a relation that is too flat, while a model with an environmentally-varying CFE but without an upper mass function truncation predicts median GC metallicities that are too high. Of the four formation models, only the fiducial E-MOSAICS model (with both an environmentally-varying CFE and |$M_{\rm c,\ast }$|⁠) reasonably matches the slope and normalization observed relation, except for the excess of high metallicity galaxies in galaxies with |$M_\ast \sim 10^{10.5} \, \rm {M}_{\odot }$| (which appears in all models).

• We use Gaussian mixture modelling (Muratov & Gnedin 2010) to test the predicted fraction of galaxies with bimodal distributions (Section 3.4). We find that 37 ± 2 per cent of simulated GC distributions are bimodal, compared with |$44^{+10}_{-9}$| per cent for observed galaxies. The mean metallicities of the metal-poor and metal-rich peaks, as well as the fraction of metal-rich GCs (for |$M_\ast \gt 10^{9.5} \, \rm {M}_{\odot }$|⁠), are also in good agreement with the observed relations.

• We demonstrate that bimodal GC metallicity distributions in E-MOSAICS galaxies are created in two ways (Section 3.4): variations in cluster formation properties as a function of metallicity that form a bimodal distribution, and cluster disruption generating a bimodal distribution from an initially-unimodal distribution. For high-mass galaxies in the simulation (⁠|$M_\ast \gt 10^{10} \, \rm {M}_{\odot }$|⁠), most initial GC metallicity distributions are unimodal, and thus bimodality occurs largely as a result of cluster disruption. Based on the metallicity distribution of stars in the Milky Way which may originate from disrupted GCs (Schiavon et al. 2017), we suggest that the bimodal GC metallicity distribution of the Milky Way is similarly a result of cluster disruption, rather than formation (Section 4.1).

Overall, the fiducial E-MOSAICS model reproduces many features of observed GC metallicity distributions, with a few identified shortcomings. Addressing the issues with insufficient GC disruption and overabundance of metal-rich GCs (for galaxies with |$M_\ast \sim 10^{10.5} \, \rm {M}_{\odot }$|⁠) requires progress on two fronts: higher resolution and more detailed galaxy formation models (particularly with regard to modelling the cold ISM; e.g. Reina-Campos et al. 2022b); and more realistic subgrid models for the dynamical evolution of GCs (currently the size evolution of GCs is not modelled, which is important for understanding how GCs respond to tidal shocks, e.g. Gieles et al. 2006a; Gieles & Renaud 2016; Webb, Reina-Campos & Kruijssen 2019; Martinez-Medina et al. 2022).

We have shown that a scenario in which old GCs and young star clusters form and evolve via the same mechanisms is able to explain the metallicity distributions observed for GCs in nearby galaxies. In view of this result, we conclude that different formation mechanisms or truncated formation epochs are not required to separately explain metal-poor and metal-rich GCs (Section 4.2). Instead, GCs can be considered as the evolved, surviving analogues of young star clusters observed forming in galaxies at the present day.

ACKNOWLEDGEMENTS

We thank Chris Usher for helpful discussions and providing SLUGGS survey data, Katja Fahrion for providing Fornax cluster data and the referee, Prof Oleg Gnedin, for helpful comments. JP is supported by the Australian government through the Australian Research Council’s Discovery Projects funding scheme (DP200102574). JMDK gratefully acknowledges funding from the German Research Foundation (DFG - Emmy Noether Research Group KR4801/1-1). JMDK and STG gratefully acknowledge funding from the European Research Council (ERC-StG-714907, MUSTANG). COOL Research DAO is a Decentralised Autonomous Organisation supporting research in astrophysics aimed at uncovering our cosmic origins. RAC is supported by the Royal Society. This work used the DiRAC Data Centric system at Durham University, operated by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment was funded by BIS National E-infrastructure capital grant ST/K00042X/1, STFC capital grants ST/H008519/1 and ST/K00087X/1, STFC DiRAC Operations grant ST/K003267/1 and Durham University. DiRAC is part of the National E-Infrastructure. The work also made use of high performance computing facilities at Liverpool John Moores University, partly funded by the Royal Society and LJMU’s Faculty of Engineering and Technology.

DATA AVAILABILITY

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

Footnotes

REFERENCES

APPENDIX: MEDIAN METALLICITIES WITH RADIUS LIMITS

In Section 3.3 we compare the median GC population metallicities of galaxies for the four E-MOSAICS GC formation models (Fig. 7). In Section 3.1 we also show that observational aperture limits for Fornax cluster galaxies have a significant effect on the predicted GC metallicity distributions for massive galaxies (Figs 2 and 3). Therefore, in Fig. A1 we repeat the comparison of median GC metallicities with a |$20 \, \rm {kpc}$| radius limit (as in Fig. 3) to demonstrate the effect of the observational aperture on the medial metallicity. In this case, the predicted median GC metallicities for massive galaxies (⁠|$M_\ast \gtrsim 5 \times 10^{10} \, \rm {M}_{\odot }$|⁠) are in better agreement with Fornax/Virgo clusters and SLUGGS survey galaxies.