Young star cluster populations in the E-MOSAICS simulations

TL;DR

This study uses the E-MOSAICS simulations to analyze young star cluster populations, demonstrating that environmentally-varying formation models align well with observed data across different redshifts and galaxy environments.

Contribution

It introduces and validates a model of star cluster formation that varies with environment, improving the understanding of cluster formation processes in galaxies.

Findings

Models accurately reproduce observed cluster populations in nearby galaxies.

High-redshift clusters have larger maximum masses due to higher gas fractions.

Environmentally-varying models outperform constant-formation-fraction models in matching observations.

Abstract

We present an analysis of young star clusters (YSCs) that form in the E-MOSAICS cosmological, hydrodynamical simulations of galaxies and their star cluster populations. Through comparisons with observed YSC populations, this work aims to test models for YSC formation and obtain an insight into the formation processes at work in part of the local galaxy population. We find that the models used in E-MOSAICS for the cluster formation efficiency and high-mass truncation of the initial cluster mass function ($M_\mathrm{c,\ast}$) both quantitatively reproduce the observed values of cluster populations in nearby galaxies. At higher redshifts ($z \geq 2$, near the peak of globular cluster formation) we find that, at a constant star formation rate (SFR) surface density, $M_\mathrm{c,\ast}$ is larger than at $z=0$ by a factor of four due to the higher gas fractions in the simulated high-redshift galaxies. Similar processes should be at work in local galaxies, offering a new way to test the models. We find that cluster age distributions may be sensitive to variations in the cluster formation rate (but not SFR) with time, which may significantly affect their use in tests of cluster mass loss. By comparing simulations with different implementations of cluster formation physics, we find that (even partially) environmentally-independent cluster formation is inconsistent with the brightest cluster-SFR and specific luminosity-$\Sigma_\mathrm{SFR}$ relations, whereas these observables are reproduced by the fiducial, environmentally-varying model. This shows that models in which a constant fraction of stars form in clusters are inconsistent with observations.

Figures

Equations

$$M_{\mathrm{V}}^{\mathrm{brightest}}=-1.91(\pm 0.09)\times \log \frac{\mathrm{SFR}}{{\mathrm{M}}_{\text{sun}}{\mathrm{yr}}^{-1}}-12.58(\pm 0.13),$$

Full text

††pubyear: 2019††pagerange: Young star cluster populations in the E-MOSAICS simulations–C

1 Introduction

Star clusters are a natural by-product of the star formation process (for recent reviews, see Longmore et al., 2014; Kruijssen, 2014; Adamo & Bastian, 2018; Krumholz et al., 2018). Young star clusters (YSCs) are observed in all star-forming galaxies for which they can be resolved (e.g. Larsen & Richtler, 1999); with the resolving power of the Hubble Space Telescope they can be detected out to distances of $\sim$100 Mpc (e.g. Adamo et al., 2010; Fedotov et al., 2011). This observability makes YSCs important tracers of the star formation process in galaxies. The most massive YSCs are also thought to be young analogues of globular clusters (GCs) (Portegies Zwart, McMillan & Gieles, 2010; Kruijssen, 2015; Forbes et al.,, 2018), therefore understanding the formation of YSCs may reveal important clues about the formation of GCs.

In recent years, observational studies have established the close connection between the properties of YSC populations and the intensity of star formation in their host galaxies (for a recent review, see Adamo & Bastian, 2018). The fraction of stars formed in bound clusters, i.e. the cluster formation efficiency (CFE or $\Gamma$, Bastian, 2008), correlates with the star formation rate (SFR) surface density $\Sigma_{\mathrm{SFR}}$ (Goddard et al., 2010; Adamo et al., 2011; Silva-Villa et al., 2013; Adamo et al., 2015; Johnson et al., 2016). At the low-mass end of their mass range, YSCs are observed to have a power-law mass function with an exponent $\beta\approx-2$ (Zhang & Fall, 1999; Bik et al., 2003; Gieles et al., 2006b; McCrady & Graham, 2007; Dowell et al., 2008; Fall & Chandar, 2012; Baumgardt et al., 2013). Both of these observations are consistent with being a natural outcome of star formation in a hierarchical gas distribution, with clusters forming in the densest regions of the gas (Elmegreen & Efremov, 1997; Efremov & Elmegreen, 1998; Elmegreen & Elmegreen, 2001; Kruijssen, 2012). However, at the high-mass end, evidence suggests that clusters form with a high-mass exponential truncation to the power-law mass function ($M_{\mathrm{c,\ast}}$) that scales with $\Sigma_{\mathrm{SFR}}$ (Gieles et al., 2006a; Gieles, 2009; Larsen, 2009; Portegies Zwart et al., 2010; Johnson et al., 2017). The observed relation between the magnitude of the brightest cluster in a population ($M_{\mathrm{V}}^{\mathrm{brightest}}$) and the SFR of the galaxy (Billett, Hunter & Elmegreen, 2002; Larsen, 2002) also implies an upper truncation mass, as it cannot be simply explained by a statistical (size-of-sample) effect with a power-law mass function (Bastian, 2008). Instead, the high-mass end of the initial cluster mass function is likely set by a combination of galactic dynamics and stellar feedback (Reina-Campos & Kruijssen, 2017).

The dependence of star cluster formation on galactic scale properties means that modelling the formation of realistic star cluster populations also requires modelling the formation and evolution of galaxies and their environment. In part for this reason, simulations of YSC populations have lagged behind the progress of observations. For computational reasons, most works modelling YSC populations focus on isolated or merging galaxies in idealised, non-cosmological simulations (e.g. Li et al., 2004; Bournaud et al., 2008; Kruijssen et al., 2011, 2012; Renaud et al., 2015; Maji et al., 2017). For the same reason, simulations in a cosmological context also largely focus on high-redshift conditions (e.g. Li et al., 2017; Kim et al., 2018), which does not allow for direct comparisons to present-day galaxies. Moreover, most studies do not investigate populations of galaxies, meaning scaling relations between YSC and galaxy properties generally cannot be compared comprehensively to observations.

In this work we investigate the YSC populations in simulations from the MOdelling Star cluster population Assembly In Cosmological Simulations within EAGLE (E-MOSAICS) project (Pfeffer et al., 2018; Kruijssen et al., 2019a). E-MOSAICS couples the MOSAICS model for star cluster formation and evolution (Kruijssen et al., 2011; Pfeffer et al., 2018) to the Evolution and Assembly of GaLaxies and their Environments (EAGLE) galaxy formation model (Crain et al., 2015; Schaye et al., 2015), therefore capturing both the evolution of the galaxies and their environment, as well as the formation and evolution of their star cluster populations. The E-MOSAICS project aims to test the origin and evolution of GC populations within a YSC-based cluster formation scenario (Pfeffer et al., 2018, 2019; Reina-Campos et al., 2018, 2019; Usher et al., 2018) and the use of star clusters as tracers of galaxy formation and assembly (Kruijssen et al., 2019a, b; Hughes et al., 2019). In the fiducial cluster formation model, star cluster populations are fully described by the local, environmentally-varying CFE (Kruijssen, 2012) and cluster truncation mass (Reina-Campos & Kruijssen, 2017). Though the analytical formulations of both models have previously been tested against observations (Kruijssen, 2012; Silva-Villa et al., 2013; Adamo et al., 2015; Johnson et al., 2016; Reina-Campos & Kruijssen, 2017; Messa et al., 2018b), the local formulation of the models and their coupling to hydrodynamical simulations through the MOSAICS model has not been systematically compared with observed YSC populations. The simulations allow for each component of the model to be switched off, such that their role in the formation of YSC populations, and the variance with galaxy properties, can be assessed. This work also serves as a means to validate the E-MOSAICS cluster formation model and thereby motivate its application to GC populations.

This paper is structured as follows. In Section 2, we briefly summarize the E-MOSAICS simulations and introduce a new set of $12.5$ comoving Mpc (cMpc) periodic volumes, for which this paper presents the first results. In Section 3, we present the results from the simulations and compare them to observations, for the CFE-$\Sigma_{\mathrm{SFR}}$ relation (Section 3.1), $M_{\mathrm{c,\ast}}$-$\Sigma_{\mathrm{SFR}}$ relation (Section 3.2), power-law indices of the mass functions (Section 3.3), specific luminosities (Section 3.4), $M_{\mathrm{V}}^{\mathrm{brightest}}$-SFR relations (Section 3.5) and cluster age distributions (Section 3.6). Finally, we summarize and discuss our conclusions in Section 4.

2 Methods

In this section we briefly describe the E-MOSAICS model and simulation suite, the selection of simulated galaxies and their star clusters and the methods for analysing the simulations. A full description of the MOSAICS model, the coupling of MOSAICS to the EAGLE model, along with extensive testing of the subgrid models, is given by Pfeffer et al. (2018), and the extension to the full suite of 25 zoom-in simulations is presented in Kruijssen et al. (2019a).

2.1 The E-MOSAICS simulations

The E-MOSAICS project (Pfeffer et al., 2018; Kruijssen et al., 2019a) is a suite of cosmological, hydrodynamical simulations of galaxy formation in the $\Lambda$ cold dark matter cosmogony that couples the MOSAICS model for star cluster formation and evolution (Kruijssen et al., 2011; Pfeffer et al., 2018) to the EAGLE model of galaxy formation and evolution (Schaye et al., 2015; Crain et al., 2015). The simulations are run with a highly modified version of the $N$-body, smoothed particle hydrodynamics code gadget3 (last described by Springel, 2005). Bound galaxies (subhaloes) were identified within the simulations using the subfind algorithm (Springel et al., 2001; Dolag et al., 2009), in the same manner as in the EAGLE simulations (for details see Schaye et al., 2015). EAGLE includes subgrid routines describing radiative cooling (Wiersma et al., 2009a), star formation (Schaye & Dalla Vecchia, 2008), stellar evolution and mass loss (Wiersma et al., 2009b), the seeding and growth of black holes (BHs) via gas accretion and BH-BH mergers (Rosas-Guevara et al., 2015), and feedback associated with star formation and BH growth (Booth & Schaye, 2009). As current cosmological simulations lack the resolution and physics necessary to compute the feedback efficiencies from first principles, the stellar and active galactic nuclei feedback parameters are calibrated such that the simulations of cosmologically representative volumes reproduce the galaxy stellar mass function, galaxy sizes and BH masses at $z\approx 0$. The EAGLE simulations successfully reproduce a range of galaxy properties, including the stellar masses (Furlong et al., 2015) and sizes (Furlong et al., 2017) of galaxies, their luminosities and colours (Trayford et al., 2015), their cold gas properties (Lagos et al., 2015, 2016; Bahé et al., 2016; Marasco et al., 2016; Crain et al., 2017), and the properties of circumgalactic and intergalactic absorption systems (Rahmati et al., 2015; Rahmati et al., 2016; Oppenheimer et al., 2016; Turner et al., 2016, 2017). The simulations also largely reproduce the cosmic star formation rate density and relation between specific star formation rate and galaxy mass (Furlong et al., 2015). The simulations are therefore ideal for comparisons with observed galaxy populations.

In the MOSAICS model, star clusters are treated as a subgrid component of star particles in the simulation (Kruijssen et al., 2011). Star clusters are therefore ‘attached’ to star particles, such that they adopt the properties of the host particle (i.e. positions, velocities, ages, abundances). In a newly formed star particle, a population of star clusters may be formed with properties that depend on the cluster formation model. The model describes cluster formation in terms of two parameters: the CFE ($\Gamma$) and the high-mass exponential truncation of the Schechter (1976) cluster mass function $M_{\mathrm{c,\ast}}$ (with a power-law index of $-2$ at lower masses). Clusters are drawn from the mass function between masses of $10^{2}$ and $10^{8}\,\mathrm{M}_{\sun}$, while only clusters with masses $>5\times 10^{3}\,\mathrm{M}_{\sun}$ are evolved to reduce the memory requirements of the simulations. Each stellar particle forms (statistically) a fraction of its mass in bound clusters ($\Gamma$ times the particle mass). Thus particles may host clusters more massive than the stellar mass of the particle, and the CFE and cluster mass function are only well sampled for an ensemble of star particles. However, the total cluster and field star mass is conserved on galactic scales (see Pfeffer et al. 2018 for details).

In the E-MOSAICS suite, we consider four variations of the cluster formation model to assess the importance of each component. The fiducial cluster formation model is fully environmentally dependent. The CFE is determined by the local formulation of the Kruijssen (2012) model, which varies as a function of the local natal gas pressure. The mass function truncation mass is determined by the local formulation of the Reina-Campos & Kruijssen (2017) model, where $M_{\mathrm{c,\ast}}$ is related to the local Toomre (1964) mass. The model assumes that $M_{\mathrm{c,\ast}}$ is proportional to the mass of the molecular cloud from which clusters form (Kruijssen, 2014). As the simulations do not have the necessary physics and resolution to model molecular clouds, their (sub-grid) masses are calculated by assuming the local Toomre mass sets the maximum mass of molecular clouds, which may further decrease due to the effects of stellar feedback. In the model, the truncation mass generally increases with the natal gas pressure, but decreases in regions with high Coriolis or centrifugal forces (i.e. near the centres of galaxies).

The three other cluster formation model variations then consider environmentally independent versions of the CFE and $M_{\mathrm{c,\ast}}$ (see also Reina-Campos et al., 2019): (i) a constant CFE of $\Gamma=0.1$ with a pure power-law mass function (i.e. $M_{\mathrm{c,\ast}}=\infty$; no formation physics model); (ii) an environmentally varying CFE with $M_{\mathrm{c,\ast}}=\infty$ (CFE only model); (iii) an environmentally dependent $M_{\mathrm{c,\ast}}$ with $\Gamma=0.1$ ($M_{\mathrm{c,\ast}}$* only* model).

The simulations model several channels of mass loss for star clusters, namely stellar evolution, two-body relaxation, tidal shock driven mass loss and complete disruption by dynamical friction (for details, see Kruijssen et al., 2011; Pfeffer et al., 2018). Stellar evolutionary mass loss for clusters is proportional to that of the host stellar particle, calculated in the EAGLE model (Wiersma et al., 2009b). The mass loss rate from two-body relaxation is determined by the strength of the local tidal field, which is calculated via the eigenvalues of the tidal tensor at the location of the star particle. The tidal shock mass loss is also calculated via the tidal tensors, based on the derivations of Gnedin, Hernquist & Ostriker (1999) and Prieto & Gnedin (2008). Star clusters that reach a mass below $100\,\mathrm{M}_{\sun}$ are assumed to be fully disrupted. Additionally, the removal of star clusters due to dynamical friction is treated in post-processing and applied at every snapshot (though this mechanism is mainly important for massive, old clusters and has little effect on young cluster populations).

In this work, we use the 25 zoom-in simulations focussed on Milky Way-mass haloes (Pfeffer et al., 2018; Kruijssen et al., 2019a) and a new set of E-MOSAICS simulations of a $12.5\,{\rm cMpc}$ periodic volume (L012N0376). All simulations were performed using a Planck Collaboration et al. (2014) cosmology, the ‘recalibrated’ EAGLE model (see Schaye et al., 2015) and initial baryonic particle masses of $\approx 2.25\times 10^{5}\,\mathrm{M}_{\sun}$. The 25 zoom-in simulations of Milky Way-mass haloes ($M_{\mathrm{vir}}\approx 10^{12}\,\mathrm{M}_{\sun}$) were drawn from the high resolution EAGLE simulation of a $25\,{\rm cMpc}$ volume (Recal-L025N0752) and resimulated in a zoom-in fashion with the E-MOSAICS model (see Pfeffer et al., 2018; Kruijssen et al., 2019a). Ten of these zoom-in simulations were run with all four MOSAICS model variations, while the other fifteen were run only with the fiducial model. To increase the range of galaxy types and environments (mainly for galaxies with stellar masses $<10^{10}\,\mathrm{M}_{\sun}$), we also performed (for all four model variations) new simulations of a periodic volume with side length $L=12.5\,{\rm cMpc}$, using $2\times 376^{3}$ particles (i.e. the EAGLE Recal-L012N0376 volume). The E-MOSAICS L012N0376 volume, and six example galaxies from the fiducial cluster formation model, are visualized in Fig. 1. In general, cluster formation is biased towards the centres of the galaxies ($\lesssim 10\,\mathrm{kpc}$ for Milky Way-mass galaxies), in regions where the natal gas densities of star formation are highest.

2.2 Galaxy and star cluster selection

We select galaxies from bound subhaloes, including both central and satellite galaxies in a halo, from both the periodic volume and zoom-in simulations at redshifts $z=0$, $0.5$, $1$ and $2$. For the zoom-in simulations, we only consider galaxies that fall within the high-resolution region of the simulations (those that reside in haloes with $<0.1$ per cent contamination by low-resolution particles at any snapshot).

Galaxies (and their bound particles) are then selected from the simulations for analysis in the following way.

• •

First, we select star particles within a $30\,\mathrm{kpc}$ three-dimensional radius (as for Schaye et al., 2015) from the centre of potential of the galaxy (i.e. the position of the most bound particle in the subhalo). This focuses the particle selection on the main galaxy and helps to exclude particles being stripped from merging satellites. Galaxies must have at least 20 star particles within this region, giving a minimum resolved stellar mass of $\approx 4\times 10^{6}\,\mathrm{M}_{\sun}$.

• •

Next, galaxies are limited to having at least 10 young ($<300\,\mathrm{Myr}$) star particles within a projected radius $R_{\mathrm{lim}}$, where the galaxy is projected such that the disc is face on (using the angular momentum of the star particles). This selection imposes a minimum total stellar mass of $>10^{7}\,\mathrm{M}_{\sun}$ at $z=0$ and $>4\times 10^{6}\,\mathrm{M}_{\sun}$ at $z=2$. We calculate $R_{\mathrm{lim}}$ as the minimum of $1.5R_{1/2}$ (the projected half-mass radius) and the radius containing 68 per cent of the recent ($<300\,\mathrm{Myr}$) star formation in the galaxy. These selections are made in order to approximate the typical footprints for observations of nearby star-forming galaxies (e.g. Adamo et al., 2015; Messa et al., 2018a) and to limit the projected region such that area-averaged quantities (e.g. $\Sigma_{\mathrm{SFR}}$) are focussed on the main star-forming component of each galaxy, respectively. The latter selection is important in galaxies with very centrally-concentrated star formation. Note that, because of the scale-free nature of the interstellar medium (ISM) and star formation (e.g. Elmegreen & Falgarone, 1996; Elmegreen, 2002), there is no standard definition for the star-forming area of galaxies (see also Kruijssen & Bastian, 2016, for a discussion on appropriate areas). Due to the heterogeneous coverage of observed galaxies, it is not possible to match the observational footprints directly (e.g. see fig. 1 of Larsen, 2002).

• •

Finally, we select star-forming galaxies based on their specific star formation rate (sSFR) measured within $1.5R_{1/2}$. Following Bourne et al. (2017), we use their equation 6 to define the star-forming galaxy sequence as a function of redshift, but set a constant sSFR for galaxies with $M_{\ast}\leq 10^{10}\,\mathrm{M}_{\sun}$ and apply a vertical shift to lower sSFRs (by setting $b_{0}=-10.2$ and $b_{1}=2.3$ in their equation 6) to match the EAGLE main sequence (which predicts slightly lower sSFRs than observed, see Furlong et al., 2015). We then select galaxies with sSFRs that do not fall more than $0.5$ dex below the sequence. At $z=0$ and stellar masses of $M_{\ast}\leq 10^{10}\,\mathrm{M}_{\sun}$, this selects galaxies with $\mathrm{SFR}/M_{\ast}>4\times 10^{-10}\,\mathrm{yr}^{-1}$.

Star clusters are selected in galaxies following the same criteria as for the star particles to which they are attached. With the exception of the CFE (which is calculated from the total initial mass in clusters) and when fitting initial cluster mass functions, we apply a mass limit for evolved clusters of $M>5\times 10^{3}\,\mathrm{M}_{\sun}$. Though this limit is necessary due to instantaneous disruption of low-mass clusters in the simulations, it is consistent with those imposed in observations of YSCs in nearby galaxies (depending on distance to the galaxy and the upper age limit for clusters; e.g. Annibali et al., 2011; Adamo et al., 2015; Johnson et al., 2017; Messa et al., 2018a; Cook et al., 2019).

For the $z=0$ snapshot, these selection criteria give us a sample of 153 galaxies with stellar masses between $2\times 10^{7}$ and $4\times 10^{10}\,\mathrm{M}_{\sun}$ (median $3.7\times 10^{8}\,\mathrm{M}_{\sun}$), SFRs between $8\times 10^{-3}$ and $3\,\mathrm{M}_{\sun}\,\mathrm{yr}^{-1}$ (median $0.04\,\mathrm{M}_{\sun}\,\mathrm{yr}^{-1}$) and $\Sigma_{\mathrm{SFR}}$ between $10^{-4}$ and $0.3\,\mathrm{M}_{\sun}\,\mathrm{yr}^{-1}\,\mathrm{kpc}^{-2}$ (median $2\times 10^{-3}\,\mathrm{M}_{\sun}\,\mathrm{yr}^{-1}\,\mathrm{kpc}^{-2}$). Of these, 39 galaxies have $>50$ YSCs (ages $<300\,\mathrm{Myr}$ and initial masses $>5\times 10^{3}\,\mathrm{M}_{\sun}$) within $R_{\mathrm{lim}}$.

2.3 Analysis

All cluster and SFR-related quantities (SFR, $\Sigma_{\mathrm{SFR}}$) are calculated for clusters and star particles with ages $<300\,\mathrm{Myr}$ at the time of the relevant snapshot, with the exception of Section 3.5 (which investigates the $M_{\mathrm{V}}^{\mathrm{brightest}}$-SFR relation) and Section 3.6 (which investigates cluster age distributions). For clusters, this age limit is similar to observational studies for which YSC populations are typically only complete (in mass) below ages of a few hundred megayears (depending on the mass limit). Star formation rates for the simulated galaxies are calculated directly from the mass in star particles formed over this time period. Observational studies often use SFR tracers (e.g. H$\alpha$ or UV flux, stellar counts) which are sensitive to timescales $\lesssim 100\,\mathrm{Myr}$ (Kennicutt & Evans, 2012; Haydon et al., 2018).

Projected galaxy quantities ($\Sigma_{\mathrm{SFR}}$, $\Sigma_{\mathrm{gas}}$, $\Sigma_{\ast}$) are calculated within the surface area given by the projected radius $R_{\mathrm{lim}}$ (i.e. area-weighted surface densities). This method follows most observational studies which use the same procedure, though it remains sensitive to the region over which the properties are measured (e.g. particularly if star formation is highly concentrated or substructured; see also Johnson et al. 2016, who apply a mass-weighted method, and Appendix A).

In Section 3.5 we compare YSC properties against the SFR and sSFR of the galaxy. For this comparison, we calculate all properties within $1.5R_{1/2}$ so as not to bias the sSFR measurement for cases with a very small $R_{\mathrm{lim}}$.

In Sections 3.2 and 3.3 we fit Schechter (1976) mass functions to the simulated YSC populations. We follow a similar analysis to that used in observational studies (e.g. Johnson et al., 2017; Messa et al., 2018b). For each population, we use the Markov chain Monte Carlo (MCMC) code pymc (Fonnesbeck et al., 2015) to sample the posterior probability distribution of the Schechter truncation mass. For each population we sample the truncation mass in $\log$-space with a uniform prior between $5\times 10^{3}\,\mathrm{M}_{\sun}$ and $10^{8}\,\mathrm{M}_{\sun}$ (the lowest mass YSC we consider and $\sim$30 times the mass of the most massive YSC at $z=0$ in our study, respectively). When fitting for the initial cluster masses, we fix the power-law index of the mass function to $\beta=-2$, the input index in the cluster models. When fitting for the final (evolved) masses of clusters, we allow the power-law index to vary, and sample the index with a uniform prior between $-3$ and $-0.5$. For each population, we perform 10,000 iterations with 1,000 burn-in steps.

For Sections 3.4 and 3.5, Johnson $U$ and $V$-band luminosities were calculated for clusters assuming simple stellar populations using the clusters’ age, metallicity and mass in combination with the Flexible Stellar Population Synthesis (FSPS) model (Conroy, Gunn & White, 2009; Conroy & Gunn, 2010), 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 and assuming the default FSPS parameters. Mass-to-light ratios for the clusters were calculated by linearly interpolating the luminosities and relative stellar masses from the grid in ages and total metallicities $\log(Z/\,\mathrm{Z}_{\sun})$. Note that we do not include extinction in these estimates, as most observational studies correct for this effect.

3 Results

3.1 Cluster formation efficiency

In Fig. 2, we first test the CFE-$\Sigma_{\mathrm{SFR}}$ relation in the fiducial cluster formation model for stars and clusters younger than $300\,\mathrm{Myr}$ at $z=0$. Note that this is not strictly a prediction of the simulations, since the Kruijssen (2012) model, for which the galaxy-scale version has previously been tested against observations, was adopted for the CFE model in the simulations. However, the section serves as a validation of the implementation in terms of local variables within the E-MOSAICS model and its extension to galaxy-wide scales within the simulations, and enables the testing of effects that may induce scatter in the relation. To calculate the CFE in the simulations, we sum the total mass of clusters formed111For a small number of particles where $M_{\mathrm{c,\ast}}<100\,\mathrm{M}_{\sun}$ (less than the lower mass for the cluster mass function) we assume that no clusters were formed. in the star particles in the region of interest (i.e. the CFE at formation which includes any stochasticity in sampling initial cluster masses from the mass function, not the values calculated at the particle level from the natal gas pressure). This value does not include the effect of cluster mass loss (which will lower the measured value of the CFE) or any observational uncertainties associated with cluster detection and measurement of the SFR. For the simulated galaxies, $\Sigma_{\mathrm{SFR}}$ is also measured over the same $300\,\mathrm{Myr}$ timescale as the CFE. Our results are not sensitive to the exact age limit used, as we find consistent results222Though with significantly fewer galaxies in the $1$-$10\,\mathrm{Myr}$ age range, since this range is generally poorly resolved due to the resolution limits of the simulations. when using ranges of $1$-$10\,\mathrm{Myr}$, $1$-$100\,\mathrm{Myr}$ and $10$-$100\,\mathrm{Myr}$.

Fig. 2 shows the galaxy-averaged CFE for all galaxies in the simulations with the fiducial cluster formation model. For comparison, in the figure we also include observed CFEs from Goddard et al. (2010), Adamo et al. (2011), Annibali et al. (2011), Silva-Villa & Larsen (2011), Cook et al. (2012, using the results for ages $<100\,\mathrm{Myr}$ and excluding galaxies with only upper limits for the CFE), Ryon et al. (2014), Adamo et al. (2015), Hollyhead et al. (2016), Johnson et al. (2016), Ginsburg & Kruijssen (2018) and Messa et al. (2018b). Note that the CFE for NGC 4449 is a lower limit (Annibali et al., 2011). The simulations show a similar level of scatter in the CFE at a given $\Sigma_{\mathrm{SFR}}$ ($\sim 0.25$ dex) to measurements of observed galaxies. In the figure there is a galaxy mass gradient along the CFE-$\Sigma_{\mathrm{SFR}}$ relation, such that the more massive galaxies generally have a higher $\Sigma_{\mathrm{SFR}}$ and CFE. This result is expected as, assuming pressure equilibrium in the galaxies, larger galaxy masses (therefore deeper potentials) result in higher characteristic ISM pressures, and thus higher $\Sigma_{\mathrm{SFR}}$ and CFE. However, it is important to note that the volumes of the simulations are not large enough to capture rare objects, such as rapidly star-forming dwarf galaxies with high star-formation densities and CFEs (e.g. blue compact dwarfs, Adamo et al., 2011). Additionally, the selection in sSFR for star-forming galaxies (Section 2.2) selects against high-mass galaxies with low star-formation densities and CFEs.

For the projected version of the Kruijssen (2012) model (dashed lines in the figure), the gas surface density ($\Sigma_{\mathrm{gas}}$, which is the fundamental quantity in the model setting the CFE, see e.g. Kruijssen & Bastian 2016 and Ginsburg & Kruijssen 2018) is converted to a SFR surface density assuming the Kennicutt–Schmidt star formation relation (Kennicutt, 1998)333Note that for consistency with the simulations, we adopt $\Sigma_{\mathrm{SFR}}=1.515\times 10^{-4}\,\mathrm{M}_{\sun}\,\mathrm{yr}^{-1}\,\mathrm{kpc}^{-2}\,(\Sigma_{\mathrm{gas}}/1\,\mathrm{M}_{\sun}\,\mathrm{kpc}^{-2})^{1.4}$, consistent with a Chabrier (2003) initial stellar mass function (see Schaye et al., 2015).. Additionally, star formation in EAGLE is implemented following the Kennicutt–Schmidt relation, rewritten as a pressure law (Schaye et al., 2015). Therefore, the naive expectation is that the simulations should broadly reproduce the (shifted) Kruijssen (2012) relation in Fig. 2 (black dashed line), where the relation has been shifted to higher $\Sigma_{\mathrm{SFR}}$ by $\approx 0.6$ dex to account for the change from the Krumholz & McKee (2005) $P$-$\Sigma_{\mathrm{SFR}}$ relation to the (Schaye & Dalla Vecchia, 2008) relation used in EAGLE. At $\Sigma_{\mathrm{SFR}}\gtrsim 5\times 10^{-3}\,\mathrm{M}_{\sun}\,\mathrm{yr}^{-1}\,\mathrm{kpc}^{-2}$, the simulated galaxies broadly follow the expected relation (black dashed line), but are generally shifted to slightly lower $\Sigma_{\mathrm{SFR}}$, with most galaxies falling between the fiducial (grey dashed line) and shifted (black dashed line) Kruijssen (2012) relations. For the modelling in this work, the fundamental relation is between the CFE and natal gas pressure, and therefore some amount of uncertainty in the CFE-$\Sigma_{\mathrm{SFR}}$ relation arises simply due to the uncertainty in the $P$-$\Sigma_{\mathrm{SFR}}$ relation. We discuss this point in further detail in Appendix A, where we show that the offset from the expected CFE-$\Sigma_{\mathrm{SFR}}$ relation is due to galaxies being offset from the expected $P$-$\Sigma_{\mathrm{SFR}}$ (see also the discussion below).

At lower surface densities ($\Sigma_{\mathrm{SFR}}\lesssim 5\times 10^{-3}\,\mathrm{M}_{\sun}\,\mathrm{yr}^{-1}\,\mathrm{kpc}^{-2}$), the simulations show a higher CFE at a given $\Sigma_{\mathrm{SFR}}$ than the Kruijssen (2012) relation, for which the cause may be multifold. Firstly, this can be caused by highly concentrated or substructured star formation, such that star and cluster formation largely occurs in a much smaller area compared to the area for which $\Sigma_{\mathrm{SFR}}$ is calculated, which will lower the measurement for $\Sigma_{\mathrm{SFR}}$ at a given CFE. Secondly, at low $\Sigma_{\mathrm{SFR}}$ there is a physical effect that increases the CFE at lower metallicities (i.e. in lower mass galaxies) due to the metallicity-dependent density threshold for star formation implemented in EAGLE444Note that this effect of increasing CFE with metallicity is not expected to occur at high $\Sigma_{\mathrm{SFR}}$, since the densities of star-forming gas in this regime are well above the metallicity-dependent threshold.. This threshold is included to model the effect of the thermogravitational collapse of warm, photoionized interstellar gas into a cold, dense phase, which is expected to occur at lower densities and pressures in metal-rich gas (Schaye, 2004). This higher density threshold at lower metallicities results (through the lower density limit for star formation imposed by the polytropic equation-of-state implemented at high gas densities) in higher pressures of star formation at a given $\Sigma_{\mathrm{SFR}}$, and therefore in higher CFEs (see fig. 3 in Pfeffer et al., 2018). Finally, variations in the natal pressure-$\Sigma_{\mathrm{SFR}}$ relation in the galaxies will, in turn, lead to variations in the CFE-$\Sigma_{\mathrm{SFR}}$ relation, through the dependence of the CFE on natal gas pressure. Such variations can be driven by random fluctuations within the galaxies (which may be most important at low $\Sigma_{\mathrm{SFR}}$), or physical variations due to differing contributions of the gravity of stars to the mid-plane gas pressure in galaxies (i.e. variations in $\phi_{P}$, see also Appendix C). We test the importance of these effects in Appendix A, finding the dominant effect to be the use of too large an area in the calculation for $\Sigma_{\mathrm{SFR}}$ (i.e. $\Sigma_{\mathrm{SFR}}$ is systematically underestimated). This effect may be mitigated by calculating a mass-weighted surface density (see Johnson et al., 2016, and Appendix A), or by judicious aperture choice, focussing on the main region of star formation. Since most studies apply the standard area-weighted calculations, we focus on that method in this paper.

The variation of the CFE at a given $\Sigma_{\mathrm{SFR}}$ can be further investigated by comparing the radial CFE distributions within the galaxies. In Fig. 3 we show the radial CFE distributions in $2\,\mathrm{kpc}$ annuli in 23 of the 25 $L^{\ast}$ galaxies (Milky Way-mass haloes; MW16 and MW22 are quenched and do not have young clusters, thus are excluded from the figure) from the zoom-in simulations (Kruijssen et al., 2019a). For this figure, we have not applied the limit on sSFR for the galaxies (Section 2.2) in order to sample a wide range of galactic environments. The majority of measurements in the radial distributions fall along the Kruijssen (2012) relation (as expected), and galaxies with higher SFRs generally show higher CFE and $\Sigma_{\mathrm{SFR}}$. However at low $\Sigma_{\mathrm{SFR}}$ ($\lesssim 10^{-3}\,\mathrm{M}_{\sun}\,\mathrm{yr}^{-1}\,\mathrm{kpc}^{-2}$), the simulations again show significant scatter from the Kruijssen (2012) relation, with most points falling between the fiducial relation and the reinterpreted relation from Johnson et al. (2016, which uses the Kruijssen model, but assumes a $\Sigma_{\mathrm{gas}}$-$\Sigma_{\mathrm{SFR}}$ relation based on the observations of , rather than the Kennicutt-Schmidt relation). This deviation can be attributed to $\Sigma_{\mathrm{SFR}}$ being averaged over a larger area than for which star and cluster formation is occurring and variations in the natal pressure-$\Sigma_{\mathrm{SFR}}$ relation, since the natal pressure is approximately constant at a given CFE (see also Appendix A). Similarly, $\Sigma_{\mathrm{SFR}}$ for the innermost radial bin in MW13 (at $\Gamma\approx 0.5$) deviates significantly from both the ‘global’ value (square symbol) due to very central star formation that dominates the cluster formation in the galaxy (for this galaxy $R_{\mathrm{lim}}=0.75\,\mathrm{kpc}$). One galaxy, MW05, has a CFE that is significantly below other galaxies at $\Sigma_{\mathrm{SFR}}\approx 2\times 10^{-3}\,\mathrm{M}_{\sun}\,\mathrm{yr}^{-1}\,\mathrm{kpc}^{-2}$. The galaxy has a very low median cluster truncation mass at $z=0$ of $M_{\mathrm{c,\ast}}\approx 100\,\mathrm{M}_{\sun}$, meaning many star particles with $M_{\mathrm{c,\ast}}<100\,\mathrm{M}_{\sun}$ form no clusters. This is caused by the low natal gas pressure for star formation ($P/k<10^{4}\,{\rm K}\,{\rm cm}^{-3}$) and a high stellar density (high $\phi_{P}$, see the discussion in Section 3.2) in the galaxy at that epoch.

3.2 Mass function truncation

In this section we test the model for the upper exponential truncation of the cluster mass function ($M_{\mathrm{c,\ast}}$) implemented in the MOSAICS models (Reina-Campos & Kruijssen, 2017). As for the CFE, this is not strictly an independent prediction since the galaxy-scale version of the model has previously been tested against observations (Reina-Campos & Kruijssen, 2017; Messa et al., 2018b), but it serves as a test and validation of the implementation in terms of local variables within the E-MOSAICS model. However, we also provide predictions at high redshifts which may be tested with future observations.

Using the fitting procedure described in Section 2.3, for each galaxy with $>50$ clusters we fit a Schechter (1976) mass function with an upper exponential truncation mass $M_{\mathrm{c,\ast}}$ to the initial masses of clusters younger than $300\,\mathrm{Myr}$, using a fixed low-mass index of $\beta=-2$ (i.e. the input value used in the simulations). In Appendix B, we compare the resulting $M_{\mathrm{c,\ast}}$ for fitting initial cluster mass functions (with a fixed power-law index) and final (evolved) cluster mass functions (with a variable power-law index). We find that both methods generally give consistent measurements for $M_{\mathrm{c,\ast}}$, with potentially a small offset to higher initial $M_{\mathrm{c,\ast}}$ due to stellar evolutionary mass loss (a factor of $\sim 0.1$ dex).

Following the fit, we exclude galaxies for which $M_{\mathrm{c,\ast}}$ is larger than the most massive cluster in the population. In such cases, $M_{\mathrm{c,\ast}}$ is poorly constrained since cluster formation does not fully sample up to the truncation mass. For galaxies with $M_{\mathrm{c,\ast}}<\max(M)$, fits typically have $1\sigma$ uncertainties of $<0.5$ dex. For $M_{\mathrm{c,\ast}}>\max(M)$ uncertainties reach up to $\sim 2$ dex, even for populations with $>100$ clusters. Due to the initial cluster mass limit ($5\times 10^{3}\,\mathrm{M}_{\sun}$), truncation masses are typically only able to be fit above masses of a few times $10^{4}\,\mathrm{M}_{\sun}$. This limit also biases the results to galaxies ($M_{\ast}\gtrsim 10^{9}\,\mathrm{M}_{\sun}$ at $z=0$) which have a large enough population of YSCs above the mass limit to fit a mass function. At $z=\{0,0.5,1,2\}$ this leaves us with a sample of $\{33,51,65,60\}$ galaxies.

In Fig. 4 we compare the fitted $M_{\mathrm{c,\ast}}$ for the simulated galaxies at $z=0$ with results from observed nearby galaxies (described in caption). The predicted $M_{\mathrm{c,\ast}}$ for the simulated galaxies are in good agreement with the observed galaxies, falling about the relation described by the observations (Johnson et al., 2017) over the same range in $\Sigma_{\mathrm{SFR}}$. More observations are needed to test whether the scatter found in $M_{\mathrm{c,\ast}}$ for the simulated YSCs is consistent with observed YSC populations, which is possible with (e.g.) the LEGUS survey (Calzetti et al., 2015) and the PHANGS-HST survey (Lee et al., in prep.).

Overall, the Reina-Campos & Kruijssen (2017) model for $M_{\mathrm{c,\ast}}$, and its implementation in terms of local gas and dynamical properties in the E-MOSAICS cluster formation model, performs well in reproducing YSC populations in realistic galaxy formation simulations at $z=0$. Therefore we can be confident in extending the formation model to more extreme environments, such as to the epochs of GC formation. In Fig. 5 we perform the same comparison of $M_{\mathrm{c,\ast}}$ against $\Sigma_{\mathrm{SFR}}$ for simulated galaxies at redshifts of $z=\{0,0.5,1,2\}$. For reference, we also show the measurements from observed galaxies at $z=0$ in each panel. At a given $\Sigma_{\mathrm{SFR}}$ or $M_{\mathrm{c,\ast}}$, the typical stellar mass of galaxies forming clusters decreases with increasing redshift, implying that galaxies of a given stellar mass (at that epoch) can form higher mass clusters in the early universe compared to $z=0$.

At $z=0.5$, the $M_{\mathrm{c,\ast}}$-$\Sigma_{\mathrm{SFR}}$ distribution is similar to $z=0$. Due to the larger sample of galaxies at this snapshot, the distribution extends to higher $M_{\mathrm{c,\ast}}$ and $\Sigma_{\mathrm{SFR}}$, comparable with that found for the Antennae galaxies ($M_{\mathrm{c,\ast}}\sim 2\times 10^{6}\,\mathrm{M}_{\sun}$, Zhang & Fall, 1999; Jordán et al., 2007). In fact, the simulated galaxy at $z=0.5$ closest to the Antennae measurement is one of a pair of galaxies about to undergo a major merger (with stellar masses $\sim 10^{10}\,\mathrm{M}_{\sun}$), which are separated by $\approx 10\,\mathrm{kpc}$ at the time of the snapshot. Due to the infrequency of the snapshots (we output 29 between $z=20$ and $z=0$), catching a galaxy merger during its peak is extremely unlikely. At higher redshifts of $z=\{1,2\}$, near the peak of GC formation for $L^{\ast}$ galaxies in the E-MOSAICS model ($z\sim 1$-$4$, Reina-Campos et al., 2019), we find in the simulated galaxies that $M_{\mathrm{c,\ast}}$ is higher at a given $\Sigma_{\mathrm{SFR}}$ than the relation at $z=0$. Therefore, to reach a given $M_{\mathrm{c,\ast}}$, galaxies require a lower $\Sigma_{\mathrm{SFR}}$ in the early universe compared to today (by $\sim$0.5 dex at $z=2$).

This increase in $M_{\mathrm{c,\ast}}$ is caused by two effects. At late times, a higher contribution of the gravity of stars to the mid-plane gas pressure (i.e. $\phi_{P}$, Elmegreen 1989; Krumholz & McKee 2005; equations 7 and 8 in Pfeffer et al. 2018) results in a lower gas surface density (and therefore Toomre mass) at a given pressure. Additionally, the density threshold for star formation increases with decreasing metallicity in the EAGLE model (which is mainly has an effect at $\Sigma_{\mathrm{SFR}}<10^{-2}\,\mathrm{M}_{\sun}\,\mathrm{yr}^{-1}\,\mathrm{kpc}^{-2}$, see the discussion in Section 3.1), thus resulting in a higher $M_{\mathrm{c,\ast}}$. We further discuss the effect of $\phi_{P}$ in Appendix C, finding that the $\phi_{P}$ increases from $\phi_{P}\approx 1$ at $z=2$ to $\phi_{P}\approx 2.5$ at $z=0$. In Fig. 5, we show the effect of decreasing $\phi_{P}$ at higher redshifts by rescaling the fit to observed local galaxies at $z=0$ assuming the median $\phi_{P}$ from the simulations at each redshift (dashed line). The simulated galaxies agree well with the rescaled relations at each redshift, demonstrating the effect of $\phi_{P}$ on $M_{\mathrm{c,\ast}}$. Note that as $\phi_{P}\approx 1$ at $z=2$ (right panel in Fig. 5), galaxies at $z>2$ should simply follow the $z=2$ relation, since $\phi_{P}$ cannot be less than unity. This result could be tested in galaxies from the local Universe by comparing cluster formation in regions of high stellar density (high $\phi_{P}$) and low stellar density (low $\phi_{P}$) at similar $\Sigma_{\mathrm{SFR}}$.

In Fig. 5 (particularly evident at $z=\{0.5,1\}$), a number of galaxies fall well below the present-day $M_{\mathrm{c,\ast}}$-$\Sigma_{\mathrm{SFR}}$ relation, approaching the value found for the CMZ in the Milky Way (Trujillo-Gomez et al., 2019). This is driven by two effects, both due to star formation in regions of high stellar density (i.e. very central star formation in the galaxy). Firstly, such galaxies have a decreased $M_{\mathrm{c,\ast}}$ due to the higher contribution of the gravity of stars to the mid-plane gas pressure (high $\phi_{P}$, as discussed above; see Fig. 16). However, a high stellar density (high $\phi_{P}$) alone is not sufficient to account for the large decrease in $M_{\mathrm{c,\ast}}$. For example, at $z=0$ a number of galaxies are significantly elevated in $\phi_{P}$ ($\phi_{P}>3$, Fig. 16), but fall along the present-day $M_{\mathrm{c,\ast}}$-$\Sigma_{\mathrm{SFR}}$ relation (Fig. 4). The main contributing factor is due to the high Coriolis or centrifugal forces in the region of star formation555A high $\phi_{P}$ does not imply high Coriolis/centrifugal forces, but high Coriolis/centrifugal forces generally occur in regions with high $\phi_{P}$., resulting in decreased Toomre masses, and therefore $M_{\mathrm{c,\ast}}$. This is captured by the Reina-Campos & Kruijssen (2017) model in terms of a dependence on the epicyclic frequency, which is not accounted for in the simple scaling with $\Sigma_{\rm SFR}$ from Johnson et al. (2017). Though this effect is most evident at $z=\{0.5,1\}$ in the simulations, it may occur at any redshift and simply results from very central star formation. Such an effect has been observed at the centre of local Universe galaxies (Reina-Campos & Kruijssen, 2017; Messa et al., 2018b; Trujillo-Gomez et al., 2019).

3.3 Mass function index

Cook et al. (2016) found that galaxies with higher star formation rate surface densities tend to have flatter cluster luminosity function indices. They suggest that this might be a result of the cluster formation process, with higher star formation efficiencies resulting in proportionally more massive star-forming regions. In this section, we investigate how other effects, namely increased relative mass loss towards low cluster masses or the degeneracy between the power-law index and $M_{\mathrm{c,\ast}}$, may instead cause this effect. Following the method described in Section 2.3, we fit Schechter and power-law mass functions to the final (evolved) cluster populations, using a variable mass function index with a uniform prior of $-3<\beta<-0.5$ (similar to Johnson et al., 2017; Messa et al., 2018b).

In Fig. 6 we compare the cluster mass function index with the star formation rate density of the galaxy, for galaxies in the $z=\{0,0.5,1,2\}$ snapshots (in order to increase the galaxy sample and extend the range in $\Sigma_{\mathrm{SFR}}$). The upper panel shows the results for Schechter function fits, while the lower panel shows the results for pure power-law fits. The figure includes all galaxies with $>50$ clusters with evolved masses $>5\times 10^{3}\,\mathrm{M}_{\sun}$ (regardless of how well $M_{\mathrm{c,\ast}}$ is constrained in the case of Schechter fits), and therefore includes the effect of the degeneracy between $\beta$ and $M_{\mathrm{c,\ast}}$ in the fits. We find that the cluster mass function indices are flatter at higher star formation rate densities, similar to the effect found by Cook et al. (2016). This result is true for both Schechter and power-law fits, though the latter tend to find steeper mass function indices. In the simulations this is caused by two effects. At high star formation rate densities ($\Sigma_{\mathrm{SFR}}\gtrsim 10^{-1.5}\,\mathrm{M}_{\sun}\,\mathrm{yr}^{-1}\,\mathrm{kpc}^{-2}$), the mass function indices become flatter due to increased relative mass loss towards low cluster masses, primarily due to tidal shocks from dense gas. At low star formation rate densities ($\Sigma_{\mathrm{SFR}}\lesssim 10^{-1.5}\,\mathrm{M}_{\sun}\,\mathrm{yr}^{-1}\,\mathrm{kpc}^{-2}$), mass function indices may appear steeper due to low truncation masses and the degeneracy between the index and the truncation mass (galaxies with $\beta<-2$ generally have poorly constrained $M_{\mathrm{c,\ast}}$). This crossing point, where galaxies typically fall above or below an index of $-2$, depends on the lower cluster mass limit; higher or lower mass limits result in higher or lower crossing points in $\Sigma_{\mathrm{SFR}}$, respectively. When using a lower cluster mass limit of $10^{4}\,\mathrm{M}_{\sun}$, rather than $5\times 10^{3}\,\mathrm{M}_{\sun}$, the crossing point shifts to higher $\Sigma_{\mathrm{SFR}}$ by $\approx 0.3$ dex. Observed local galaxies at low $\Sigma_{\mathrm{SFR}}$ ($\sim 10^{-3}$-$10^{-2}\,\mathrm{M}_{\sun}\,\mathrm{yr}^{-1}\,\mathrm{kpc}^{-2}$) are also consistent with the power-law indices for the cluster mass function found in this work ($\beta\approx-2$, e.g. M31, Johnson et al. 2017; M51, Messa et al. 2018a).

The mass function indices at a given $\Sigma_{\mathrm{SFR}}$ for $z=0$ and $z>0$ galaxies are generally consistent. However, at $\Sigma_{\mathrm{SFR}}\lesssim 10^{-2}\,\mathrm{M}_{\sun}\,\mathrm{yr}^{-1}\,\mathrm{kpc}^{-2}$, galaxies at $z=0$ tend to have steeper indices due to their smaller $M_{\mathrm{c,\ast}}$ (Section 3.2).

Since their methods differ from ours (they fit luminosity functions and have a variable lower cluster luminosity limit between galaxies), we cannot make a direct quantitative comparison to the results from Cook et al. (2016). Additionally, the simulations and observations are largely biased to different star formation rate densities ($\Sigma_{\mathrm{SFR}}\gtrsim 10^{-2.5}\,\mathrm{M}_{\sun}\,\mathrm{yr}^{-1}\,\mathrm{kpc}^{-2}$ and $\Sigma_{\mathrm{SFR}}\lesssim 10^{-2}\,\mathrm{M}_{\sun}\,\mathrm{yr}^{-1}\,\mathrm{kpc}^{-2}$, respectively). However, the results from the simulations suggest that the relation between cluster luminosity function index and $\Sigma_{\mathrm{SFR}}$ found by Cook et al. (2016) can be explained by the degeneracy between the mass function index and truncation $M_{\mathrm{c,\ast}}$ (at low $\Sigma_{\mathrm{SFR}}$). Similar measurements of the mass function index should be extended to higher star formation rate densities to assess and test the impact of cluster mass loss.

3.4 Specific $U$-band cluster luminosity

An empirical precursor to the CFE, the specific $U$-band cluster luminosity, $T_{L}(U)=100L_{\mathrm{clusters}}/L_{\mathrm{galaxy}}$ (i.e. the percentage of $U$-band light of a galaxy contributed by star clusters), was introduced by Larsen & Richtler (2000), who found it correlated strongly with $\Sigma_{\mathrm{SFR}}$ for observed galaxies. While the CFE is the most relevant quantity for simulations of cluster formation, it is not a direct observable and a number of caveats apply to its estimation (e.g. assumptions for and extrapolation of the cluster mass function, uncertainties in masses and ages from stellar population modelling, SFR estimations, corrections for dust, etc.). On the other hand, $T_{L}(U)$ can be directly determined from observations of galaxies (though may depend on selection criteria for YSCs) and therefore presents a useful test for models of YSC populations.

In Fig. 7 we show $T_{L}(U)$ for the simulated galaxies at $z=0$. We calculate the total luminosity of star clusters for all surviving clusters with initial masses $>5\times 10^{3}\,\mathrm{M}_{\sun}$. In order to calculate the total $U$-band luminosity of the galaxy, we assume simple stellar populations for each star particle and calculate their luminosities using the method described in Section 2.3. Both luminosities were calculated within $R_{\mathrm{lim}}$ and assuming no extinction666Adopting a basic model for extinction where all stars and clusters are embedded within an optically thick cloud until a specific age (e.g. $10\,\mathrm{Myr}$; c.f. Charlot & Fall, 2000) has no effect on the results, because extinction has the same effect on stars and clusters.. We compare the simulated galaxies in Fig. 7 with observed galaxies from Larsen & Richtler (2000) and Adamo et al. (2011). We find good agreement in the trend of $T_{L}(U)$ with $\Sigma_{\mathrm{SFR}}$ between the simulated and observed galaxies. At low $\Sigma_{\mathrm{SFR}}$ ($\sim 10^{-3}\,\mathrm{M}_{\sun}\,\mathrm{yr}^{-1}\,\mathrm{kpc}^{-2}$), $T_{L}(U)$ may be slightly underestimated in the simulations due to the instantaneous disruption of clusters with initial masses $M<5\times 10^{3}\,\mathrm{M}_{\sun}$. However, similar limitations also apply to the observed galaxies, depending on the distance to the galaxy and detectability of clusters. In the simulations, $T_{L}(U)$ is largely determined by the CFE, such that those galaxies with a high CFE also have a high $T_{L}(U)$. The scatter in $T_{L}(U)$ at fixed $\Sigma_{\mathrm{SFR}}$ (or CFE) shows no clear trends with sSFR or metallicity, and arises from temporal variations in the CFE and $\Sigma_{\mathrm{SFR}}$, as well as stochastic sampling at the high-mass end of the cluster mass function.

In Fig. 8, we quantify the effect on $T_{L}(U)$ of varying the star cluster formation physics in the simulations and show the four cluster formation models described in Section 2.1. Each model is a variation on the fiducial, environmentally-dependent model (top left panel), with environmentally-independent versions for either the truncation mass or CFE (CFE only and $M_{\mathrm{c,\ast}}$* only* models, bottom left and right panels, respectively), or both (no formation physics model, top right panel). The figure only includes galaxies from the L012N0376 volume and the first ten zoom-in simulations (MW00-MW09), i.e. simulations with all four variations of the cluster formation physics (and thus the top left panel shows fewer galaxies than in Fig. 7).

Fig. 8 clearly shows the critical role of both the CFE and $M_{\mathrm{c,\ast}}$ models in reproducing the observations of $T_{L}(U)$. With a constant CFE ($\Gamma=0.1$) and pure power-law mass function (upper right panel), the no formation physics model implies a (roughly) constant $T_{L}(U)$, and therefore cannot simultaneously reproduce galaxies of high ($>10$) and low ($<1$) specific luminosities. The CFE only model (bottom left panel) assumes an environmentally-dependent CFE and a pure power-law mass function. Due to the variation of the CFE with $\Sigma_{\mathrm{SFR}}$, the model agrees better with the observed galaxies than for the model with constant CFE. However, a variation in CFE alone (at least in the current formulation of the Kruijssen 2012 model) is also largely unable to account for galaxies with $T_{L}(U)<1$. The CFE only model predicts higher $T_{L}(U)$ than observed at $\Sigma_{\mathrm{SFR}}\sim 10^{-3}\,\mathrm{M}_{\sun}\,\mathrm{yr}^{-1}\,\mathrm{kpc}^{-2}$, but agrees with the observed galaxies at higher $\Sigma_{\mathrm{SFR}}$. The bottom right panel of Fig. 8 shows the $M_{\mathrm{c,\ast}}$* only* model, which assumes a constant $\mathrm{CFE}=0.1$ and an environmentally-dependent $M_{\mathrm{c,\ast}}$. Though the model assumes the same constant CFE as for the no formation physics model, the $M_{\mathrm{c,\ast}}$* only* model shows good agreement for galaxies with $\Sigma_{\mathrm{SFR}}\lesssim 10^{-2}\,\mathrm{M}_{\sun}\,\mathrm{yr}^{-1}\,\mathrm{kpc}^{-2}$, but underpredicts $T_{L}(U)$ at higher $\Sigma_{\mathrm{SFR}}$. In the $M_{\mathrm{c,\ast}}$* only* model, $T_{L}(U)$ is significantly lower than the assumed 10 per cent CFE due to the low $M_{\mathrm{c,\ast}}$ (at low $\Sigma_{\mathrm{SFR}}$) and the preferential formation of very low mass clusters, which are not detected777Note that $T_{L}(U)$ will therefore depend upon the lower initial cluster mass limit in the simulations (we adopt $5\times 10^{3}\,\mathrm{M}_{\sun}$). However, similar detection limits also apply for observed galaxies, depending on cluster age (e.g. Johnson et al., 2015; Adamo et al., 2015; Messa et al., 2018a).. Therefore, we conclude that environmental variations in both the CFE and $M_{\mathrm{c,\ast}}$ are necessary for reproducing the observed $T_{L}(U)$ relation.

3.5 Brightest cluster-SFR relation

3.5.1 Fiducial E-MOSAICS model

The relation between the brightest cluster and the SFR of the galaxy is an empirical relation, of which the construction does not rely on modelled quantities such as clusters ages or masses. The relation is sensitive to cluster formation physics and therefore presents a useful test of cluster formation models. Moreover, neither the star cluster nor galaxy formation physics implemented in the simulations were calibrated to reproduce the relation, thus a comparison affords an independent test of the predictions of the simulations.

In this section we consider two cases. First, we consider the relation without applying an age limit to the simulated clusters, since in a number of cases observational measurements are obtained with single-band photometry, prohibiting age measurements (e.g. Larsen, 2002). Second, we apply an upper age limit of $300\,\mathrm{Myr}$, consistent with the calculation for the SFR. For the simulations we apply a cluster luminosity limit of $M_{\mathrm{V}}<-8.2$ (assuming a metallicity of $\log(Z/\,\mathrm{Z}_{\sun})=-0.5$, similar to the metallicity of star-forming gas in a $M_{\ast}\sim 10^{8}\,\mathrm{M}_{\sun}$ galaxy in the EAGLE Recal model, Schaye et al. 2015) to reflect the $5\times 10^{3}\,\mathrm{M}_{\sun}$ initial cluster mass limit, which is similar to the luminosity limit of observational programmes (e.g. Whitmore et al., 2014). V-band mass-to-light ratios for the simulated clusters are calculated using the FSPS model (see Section 2.3). Luminosities were then determined using the current cluster masses, which include cluster mass loss.

Fig. 9 shows the predictions for the brightest cluster in the V-band as a function of galaxy SFR (both measured within $1.5R_{1/2}$) for the fiducial E-MOSAICS model, with no cluster age limit (left panel) and with an upper limit of $300\,\mathrm{Myr}$ (right panel). The best-fitting relation for the simulated galaxies from a least-squares linear regression is given in each panel and shown as the solid red line. The brightest cluster is generally consistent for both with and without an age limit. However, for some galaxies, the brightest young cluster is significantly fainter than the brightest cluster of any age, which results in a slightly flatter slope of the best-fitting relation for the $<300\,\mathrm{Myr}$ age limit. For comparison, the figure also shows the sample of observational measurements compiled by Adamo et al. (2015). Note that we have not attempted to match the sample selection for the observations, other than the selection in sSFR for the simulations, and therefore some bias between the simulated and observed galaxy populations may exist. The best-fitting relation for the observed sample of galaxies is given by

[TABLE]

shown as a dashed line. To this compilation of galaxies we also add measurements of the Milky Way and M31. For the Milky Way we assume $\mathrm{SFR}=1.9\,\mathrm{M}_{\sun}\,\mathrm{yr}^{-1}$ (Chomiuk & Povich, 2011) and the brightest cluster to be Westerlund 1, with $M_{\mathrm{V}}^{\mathrm{brightest}}\approx-11.7$ (assuming a mass of $6\times 10^{4}\,\mathrm{M}_{\sun}$ and age of 5 Myr, Mengel & Tacconi-Garman, 2007, and a V-band mass-to-light ratio from FSPS assuming a Solar metallicity). For M31 we assume $\mathrm{SFR}=0.7\,\mathrm{M}_{\sun}\,\mathrm{yr}^{-1}$ (Kang et al., 2009; Lewis et al., 2015). We take the brightest cluster in M31 (of any age) to be the globular cluster G1, with $M_{\mathrm{V}}^{\mathrm{brightest}}=-10.66$ (Galleti et al., 2004). For the brightest young cluster we use the brightest cluster from the PHAT survey, $M_{\mathrm{V}}^{\mathrm{brightest}}(<1\,\mathrm{Gyr})\approx-10.46$ (Johnson et al., 2015), using their eq. 6 to calculate a V-band magnitude and assuming a distance modulus of $24.47$ (McConnachie et al., 2005).

The slope of the best-fitting relation for the simulated galaxies ($-1.89\pm 0.15$ for no age limit, $-1.69\pm 0.2$ for clusters $<300\,\mathrm{Myr}$) is fully consistent with that of the observed galaxies ($-1.91\pm 0.09$ for the Adamo et al. 2015 sample; $-1.87\pm 0.06$ for the relation from Weidner et al. 2004). The scatter around the observed relation is also very similar for the simulations and observations, with standard deviations in $M_{\mathrm{V}}^{\mathrm{brightest}}$ from the predicted and observed relations of $\approx 0.95\,\mathrm{mag}$ and $1.01\,\mathrm{mag}$, respectively. Therefore, the observed $M_{\mathrm{V}}^{\mathrm{brightest}}$-SFR relation is reproduced by the fiducial E-MOSAICS model with an environmentally varying CFE and upper mass function truncation mass, such that star formation in environments with high natal gas pressures results in more star formation in bound clusters and up to higher cluster masses.

To investigate the origin of scatter away from the $M_{\mathrm{V}}^{\mathrm{brightest}}$-SFR relation, in Fig. 9 we colour the simulated galaxies by their sSFR. In the right panel (ages $<300\,\mathrm{Myr}$), at a fixed SFR the simulations show a gradient in sSFR, such that the galaxies with the brightest clusters also typically have the highest sSFR (or lowest galaxy masses). This trend is weaker in the left panel (no age limit) as cluster luminosities may be uncorrelated with the present SFR. We explore this further in Fig. 10, where we show the magnitude difference from the observed $M_{\mathrm{V}}^{\mathrm{brightest}}$-SFR relation (Weidner et al., 2004) compared with the sSFR for the simulated galaxies for cluster ages $<300\,\mathrm{Myr}$. Simulated galaxies with brightest clusters that are significantly fainter than the observed relation typically have the lowest sSFRs. The fitted relation crosses the zero-point (in $\Delta M_{\mathrm{V}}^{\mathrm{brightest}}$) at $\mathrm{sSFR}\approx 0.2\,\mathrm{Gyr}^{-1}$, similar to the typical sSFR (which is a weak function of galaxy mass) for galaxies on the star-forming main sequence (Schiminovich et al., 2007). By selecting galaxy populations at different constant sSFRs, a prediction of our fiducial model is that we expect to find (age-limited) $M_{\mathrm{V}}^{\mathrm{brightest}}$-SFR relations that are offset to fainter luminosities at lower sSFRs. Additionally, galaxies at higher redshifts, with higher sSFRs, should be offset to brighter $M_{\mathrm{V}}^{\mathrm{brightest}}$. The physical cause of this effect in the simulations is lower cluster mass function truncations in higher mass galaxies (which at fixed SFR have lower sSFR). These galaxies typically have lower gas mass fractions, which result in a higher $\phi_{P}$, while larger galaxy masses result in an increased importance of Coriolis/centrifugal forces in setting the Toomre mass. Both factors result in a lower cluster mass function truncation. This result could be tested with future observations in galaxies of different masses at fixed sSFR.

3.5.2 Alternative cluster formation models

The result from Fig. 9, that environmentally-dependent cluster formation reproduces the $M_{\mathrm{V}}^{\mathrm{brightest}}$-SFR relation, can be further tested by considering the alternative cluster formation model variations in the E-MOSAICS suite. In Fig. 11, we compare the $M_{\mathrm{V}}^{\mathrm{brightest}}$-SFR relation for the four cluster formation models described in Section 2.1. As in Fig. 8 (Section 3.4), each model is a variation on the fiducial, environmentally-dependent model (top left panel), with environmentally-independent versions for either the truncation mass or CFE (CFE only and $M_{\mathrm{c,\ast}}$* only* models, bottom left and right panels, respectively), or both (no formation physics model, top right panel) and only includes galaxies from simulations with all four versions of the cluster formation physics (the L012N0376 volume and first ten zoom-in simulations, MW00-MW09). Fig. 11 highlights the importance of cluster formation physics in governing the $M_{\mathrm{V}}^{\mathrm{brightest}}$-SFR relation.

In the no formation physics model (top right panel), which assumes $\Gamma=0.1$ and a $\beta=-2$ power-law mass function (i.e. with our standard cluster formation physics disabled), the simulations are inconsistent with the observed relation and recover the slope of the expected relation for a pure power-law mass function and constant CFE (dotted line in the figure; Bastian, 2008). The relation is therefore determined by a size-of-sample effect, with larger cluster populations more likely to have brighter clusters. The stochasticity at the high-mass end of the mass function induces a large scatter between galaxies at a given SFR. Also for this reason, the absolute offset in the relation is determined by the choice of CFE. However, for any choice of constant CFE, the slope of the relation will remain inconsistent with the observations.

In the CFE only model (bottom left panel), which assumes an environmentally-dependent CFE and a $\beta=-2$ power-law mass function, the simulations are again inconsistent with the observed relation. Due to the correlation of the CFE with galaxy mass (and thus SFR; see also Fig. 2), the CFE only model yields a relation that is even steeper than the no formation physics model, since low-mass galaxies with low SFRs typically have CFEs below 10 per cent. Thus at low SFRs in the figure ($\mathrm{SFR}<0.1\,\mathrm{M}_{\sun}\,\mathrm{yr}^{-1}$) many galaxies are consistent with the observed relation, while those at higher SFRs are inconsistent with observed counterparts due to the lack of a truncation in the cluster mass function. Again, there is large galaxy-to-galaxy scatter at a given SFR due to the stochasticity at the high-mass end of the mass function.

Finally, the $M_{\mathrm{c,\ast}}$* only* model, with a constant $\Gamma=0.1$ and an environmentally-dependent exponential truncation mass to the cluster mass function ($M_{\mathrm{c,\ast}}$), is shown in the bottom right panel of the figure. At low SFRs, the simulations are consistent with observed galaxies. However, since the CFE does not vary between galaxies, this model predicts a flatter slope than is observed, and thus at higher SFRs the brightest clusters are under-luminous compared to observed galaxies.

These results therefore show that an environmental dependence of both the CFE and cluster truncation mass is necessary for reproducing the observed properties of YSC populations.

3.6 Cluster age distributions

In this section we investigate the effect of time-varying star and cluster formation rates (CFRs) on cluster age distributions. The age distribution of clusters is considered to be a strong test of cluster mass-loss models, potentially enabling the discrimination between mass/environmentally dependent or independent cluster disruption (e.g. Gieles et al., 2005; Lamers et al., 2005; Whitmore et al., 2007; Bastian et al., 2009; Lamers, 2009; Chandar et al., 2010; Kruijssen et al., 2011, 2012; Bastian et al., 2012; Miholics et al., 2017, see Adamo & Bastian 2018 for a review). However, the effect of time-varying SFRs and CFRs on the interpretation of cluster age distributions has not been investigated in realistic galaxy simulations in cosmological environments.

In the upper panel of Fig. 12, we show the cluster age distributions of galaxies from the zoom-in simulations of Milky-Way mass galaxies. We apply a cluster mass limit of $>5\times 10^{3}\,\mathrm{M}_{\sun}$ at all ages and use clusters within $R_{\mathrm{lim}}$ from the centre of the potential. We limit the sample to galaxies with at least $100$ surviving clusters younger than $10^{9.25}\,\mathrm{yr}$ (similar to the number of clusters in M31 more massive than $5\times 10^{3}\,\mathrm{M}_{\sun}$, Fouesneau et al., 2014; Johnson et al., 2016), leaving 12 galaxies. For the age distributions we use eight bins with widths of 0.375 dex between $10^{6.25}$ and $10^{9.25}\,\mathrm{yr}$, where the distributions are normalised at $10^{7.5}\,\mathrm{yr}$ by fitting a power-law relation to the four bins between $10^{7}$ and $10^{8.5}\,\mathrm{yr}$. The typical timesteps for stellar particles at $z\approx 0$ are $\sim 1\,\mathrm{Myr}$, and therefore even in the youngest age bins the numerical sampling of disruption by tidal heating is adequate. We find a large scatter between the simulations, with some galaxies showing very flat age distributions and others where the number of clusters decreases rapidly at ages $>10^{8}\,\mathrm{yr}$. We discuss the cause of this scatter between galaxies below. The solid line in the figure shows the total cluster age distribution for all of the galaxies combined. For the four bins in the range $\log(\mathrm{age}/\,\mathrm{yr})=\{7,8.5\}$ (where observations are complete) we find a slope of $-0.39\pm 0.04$ for the total population, with a range between $-1.18$ and $0.17$ for individual galaxies.

We compare the cluster age distributions from the simulations with the observed distributions in M31 (Fouesneau et al., 2014; Johnson et al., 2016), M83 (Silva-Villa et al., 2014; Adamo & Bastian, 2018) and M51 (Messa et al., 2018a, b), also using a mass limit of $>5\times 10^{3}\,\mathrm{M}_{\sun}$. This sample comprises observed galaxies most similar to our sample of simulated Milky Way-mass galaxies. At ages $<10^{7}\,\mathrm{yr}$, the observations may be contaminated by unbound stellar associations (Gieles & Portegies Zwart, 2011; Bastian et al., 2012; Ward & Kruijssen, 2018), while at ages $>10^{8.5}\,\mathrm{yr}$ cluster populations are incomplete due to luminosity limits. Within the region for which observations are complete ($10^{7}$-$10^{8.5}\,\mathrm{yr}$), the age distributions of the simulations and observations are in good agreement, showing a similar level of scatter between different galaxies and consistent slopes of the age distributions ($-0.07$ and $-0.14$ for the Fouesneau et al. 2014 and Johnson et al. 2016 M31 catalogues, respectively, $-0.43$ for M51 and $-0.33$ for M83). We note that mass loss by tidal shocks is underestimated in the simulations, due to the lack of an explicit model for the cold, dense phase of the star-forming interstellar medium in the EAGLE model (see Pfeffer et al., 2018). The extent to which this affects the age distribution predictions from the simulations depends upon what cluster age this mechanism becomes important at $z\approx 0$ (which likely depends upon the local environment within which young clusters reside). However, a similar result, where cluster mass loss is mainly important after a few hundred megayears, was found in simulations of isolated galaxies by Kruijssen et al. (2011, using the same MOSAICS dynamical evolution model, but using a galaxy formation model with a simple model for the cold, dense phase of the interstellar medium, see ).

The cluster age distribution is a function of both cluster formation and evolution. Therefore, variations in the SFR or CFR will also lead to variations in the age distributions. We assess the impact of these effects in the bottom two panels of Fig. 12. In the middle panel of Fig. 12, we show the cluster age distribution normalised by the total SFR in each temporal bin for each simulated galaxy. Overall, the distributions for the galaxies show similar behaviour to the standard distributions (upper panel), with little reduction in scatter between the galaxies (the scatter about the total relation decreased from $0.43$ dex to $0.34$ dex). This is expected, since the SFR in each galaxy is relatively constant over the period investigated. Over $10^{9.2}\,\mathrm{yr}$, the SFR in each temporal bin typically varies by less than a factor of two from the median for all bins. The majority of the impact of variations in the SFR occurs at $>10^{8.5}\,\mathrm{yr}$, as can be seen by the slight reduction of scatter in the $10^{9}\,\mathrm{yr}$ bin between the upper and middle panels in the figure. The slope of the SFR-normalised age distribution is also very similar to the standard age distribution ($-0.35\pm 0.07$, with a range between $-1.14$ and $0.15$). Therefore, for the galaxies investigated, any variations in SFR induce little impact in the cluster age distributions.

In the bottom panel of Fig. 12 we show the cluster age distributions normalised by the CFR (for masses $>5\times 10^{3}\,\mathrm{M}_{\sun}$) in each temporal bin for the simulated galaxies. Again, the slope of the CFR-normalised age distribution is similar to the standard age distribution, but with a reduced range ($-0.37\pm 0.06$, with a range between $-0.75$ and $-0.09$). Unlike the SFR, normalising by the CFR has a large effect on the age distributions, significantly reducing the galaxy-to-galaxy scatter in the distributions from 0.43 dex to 0.21 dex about the total distribution. When normalised by the CFR, the age distributions follow a very similar evolution for the galaxies, with some divergence at $>10^{8}\,\mathrm{yr}$ reflecting the varying cluster mass-loss rates between the galaxies (due to, e.g., differing potentials, gas densities or galactocentric radii of cluster formation). Though the SFR varies by up to a factor of two between bins, the CFR can vary by up to a factor of six, and therefore the CFR is not simply following changes in the SFR888A similar level of variation occurs when comparing the initial age distribution, $\mathrm{d}N_{\mathrm{init}}/\mathrm{d}t$ (i.e. without including cluster mass loss), rather than $\mathrm{CFR}(M>5\times 10^{3}\,\mathrm{M}_{\sun})$, thus it is not simply due to stochasticity in the masses of clusters that form.. In addition to the SFR, the CFE (i.e. the fraction of star formation in bound clusters) may also vary with time by up to a factor of two, and thus together account for about half of the variance in the $\mathrm{CFR}(M>5\times 10^{3}\,\mathrm{M}_{\sun})$. The rest of the variation (factor of three) can be attributed to the stochasticity of cluster formation at the high-mass end of the cluster mass function (seven of the galaxies have $M_{\mathrm{c,\ast}}<10^{5}\,\mathrm{M}_{\sun}$, and therefore in these galaxies only the high-mass end of the cluster population satisfies our mass cut of $M>5\times 10^{3}\,\mathrm{M}_{\sun}$). This effect is lessened in galaxies with larger cluster populations, although even for those with $>1000$ clusters above the mass limit, the variation in the mass-limited CFR, in addition to that of the SFR and CFE, is a factor of two. For the observed galaxies in Fig. 12, M31 may be most affected by this effect since it has the fewest clusters with $M>5\times 10^{3}\,\mathrm{M}_{\sun}$, though this may be alleviated somewhat by including lower mass clusters (i.e. the completeness limit in the PHAT survey for clusters younger than $300\,\mathrm{Myr}$ is $\sim 1000\,\mathrm{M}_{\sun}$, Johnson et al., 2015).

In Fig. 12, we also show the contribution of stellar-evolutionary mass loss (dash-dotted line) to the total cluster age distribution (solid line). The stellar mass loss relation has a significantly flatter slope ($\approx-0.1$) than the total cluster age distribution. Stellar mass loss dominates for ages $\lesssim 10^{7}\,\mathrm{yr}$, with dynamical mass loss (two-body relaxation and tidal shocks) becoming dominant only at older ages.

Therefore, temporal variations in the CFE and CFR should also be considered when using cluster age distributions to test models for cluster mass loss, since they may impart non-negligible variations in the cluster age distributions.

4 Discussion and conclusions

In this paper, we present a comparison of the YSC populations in the E-MOSAICS simulations with observed populations in nearby galaxies. The aim of this work is to both test the cluster formation models in realistic simulations of galaxy formation and, by varying the cluster formation physics (the CFE and exponential truncation of the cluster mass function), obtain an insight into the formation processes at work in observed cluster populations.

We find that, due to a combination of spatially varying (non-uniform) star formation and the metallicity-dependent density threshold for star formation in the EAGLE model, the CFE in the simulated galaxies is elevated above the (input) Kruijssen (2012) relation at $\Sigma_{\mathrm{SFR}}<10^{-2.5}\,\mathrm{M}_{\sun}\,\mathrm{yr}^{-1}\,\mathrm{kpc}^{-2}$ (Section 3.1; Fig. 2). In fact, our sample of simulated galaxies could be fit by a power-law relation in CFE-$\Sigma_{\mathrm{SFR}}$ (also see Goddard et al., 2010), despite the input relation. A similar effect might be present in observed cluster populations, though this could also be caused by a non-linear $\Sigma_{\mathrm{gas}}$-$\Sigma_{\mathrm{SFR}}$ relation (Johnson et al., 2016).

For cluster populations in the simulated galaxies that can be fit with a Schechter (1976) mass function, the formulation of the Reina-Campos & Kruijssen (2017) model implemented in E-MOSAICS performs well in reproducing observed truncation masses (Section 3.2; Fig. 4). However, more observations are needed to test if the scatter in $M_{\mathrm{c,\ast}}$ seen between galaxies (i.e. at a given $\Sigma_{\mathrm{SFR}}$) in the simulations is consistent with observed galaxies. Due to the higher gas fractions of galaxies at higher redshifts, the E-MOSAICS model predicts that, at the epochs of GC formation ($z\gtrsim 1$), the $M_{\mathrm{c,\ast}}$-$\Sigma_{\mathrm{SFR}}$ relation should be elevated above the present day relation (Fig. 5).

In Section 3.3, we investigate the power-law index of the cluster mass function at low cluster masses, for evolved cluster populations. We find that in galaxies with high SFR surface densities ($\Sigma_{\mathrm{SFR}}\gtrsim 10^{-1}\,\mathrm{M}_{\sun}\,\mathrm{yr}^{-1}\,\mathrm{kpc}^{-2}$), mass function indices become flatter ($\beta>-2$) due to cluster disruption and increased relative mass loss towards low cluster masses. In galaxies with lower SFR surface densities, mass function indices are similar to the initial index ($\beta=-2$) or may potentially become steeper ($\beta<-2$) due to the degeneracy between the mass function index and upper truncation mass ($M_{\mathrm{c,\ast}}$). The results are consistent with the findings of Cook et al. (2016), who investigated observed cluster luminosity functions, suggesting that variable cluster mass function index is not necessary to explain their findings.

In Section 3.5, we test the relation between the brightest cluster within a galaxy and its SFR. This observed relation presents a useful test for the simulations since it sensitive to cluster formation physics, yet was not used to test or calibrate the cluster formation models. We find that the fiducial E-MOSAICS cluster formation model reproduces both the slope of the $M_{\mathrm{V}}^{\mathrm{brightest}}$-SFR relation and the scatter around the relation (Fig. 9). If an upper age limit ($<300\,\mathrm{Myr}$) is applied to the clusters, we find that a number of the simulated galaxies fall significantly below the observed relation, which tend to be those with low sSFRs. This result implies that the $M_{\mathrm{V}}^{\mathrm{brightest}}$-SFR relation may partially arise due to the preferential selection of galaxies on the star-forming main sequence of galaxies. Selecting samples of galaxies at fixed sSFRs may therefore result in offset $M_{\mathrm{V}}^{\mathrm{brightest}}$-SFR sequences, which is a prediction offering a new avenue to test cluster formation models.

Some previous work has suggested that young cluster formation proceeds in an environmentally-independent manner (Whitmore et al., 2007; Fall & Chandar, 2012; Chandar et al., 2015; Chandar et al., 2017). However, in Sections 3.4 and 3.5.2 we find that only the fully environmentally-dependent cluster formation (fiducial) model can explain both the $T_{L}(U)$-$\Sigma_{\mathrm{SFR}}$ and $M_{\mathrm{V}}^{\mathrm{brightest}}$-SFR relations. If cluster formation is environmentally-independent (constant CFE and a power-law mass function with $\beta=-2$), both the slope of the $M_{\mathrm{V}}^{\mathrm{brightest}}$-SFR relation and its absolute offset are inconsistent with the observed relation. Similarly, a constant CFE cannot simultaneously explain galaxies at both high and low specific $U$-band cluster luminosities [$T_{L}(U)$], while an environmental-dependence in $M_{\mathrm{c,\ast}}$ is required to explain galaxies with $T_{L}(U)<1$. We therefore conclude that an environmental dependence in cluster formation, both for the CFE and the upper truncation of the cluster mass function ($M_{\mathrm{c,\ast}}$), is required to reproduce observed young cluster populations. This shows that models in which a constant fraction of stars form in clusters are inconsistent with observations. The importance of environmentally-dependent cluster formation in reproducing GC populations has similarly been discussed in previous work with the E-MOSAICS simulations (Pfeffer et al., 2018; Usher et al., 2018; Reina-Campos et al., 2019).

Finally, we compare the cluster age distributions in 12 (out of 25) of our sample of Milky Way-mass galaxies with those in observed nearby disc galaxies (Section 3.6). We find that, for ages where observed cluster populations are generally complete (between $10^{7}$ and $10^{8.5}\,\mathrm{yr}$), the scatter in the age distributions between the simulated galaxies is similar to that in the observed galaxies. However, for the environmentally-dependent cluster formation model, the variation of the CFR (not SFR) over Myr timescales may impart significant variation in the cluster age distributions (Fig. 12). Therefore, the use of cluster ages as a discriminator of cluster mass-loss mechanisms should be approached with caution unless the variation of the CFR with time can be accounted for. The variation of the CFR is mainly driven by an underlying time-variation of the CFE. However, even for galaxies with $>1000$ clusters, additional (potentially dominant) variations may be imparted simply due to the stochasticity of cluster formation in time, with larger variations for less numerous cluster populations.

The results presented in this paper reinforce the understanding of cluster formation as an environmentally-dependent process. We make predictions for future comparisons to observations ($M_{\mathrm{c,\ast}}$ at high and low gas fractions, Section 3.2; the $M_{\mathrm{V}}^{\mathrm{brightest}}$-SFR relation at varying sSFRs, Section 3.5), which will be useful in both further testing the models implemented in E-MOSAICS, as well as cluster formation theories in general. This work also highlights the importance of realistic simulations of galaxies with a diverse range of properties and environments when testing models of cluster formation. Comparing the simulations and observations at the extremes of the galaxy population will enable a strong test of cluster formation theories, by testing them outside of the range in which they were developed, and allow further insight into the star and galaxy formation mechanisms which shape young star cluster populations across cosmic history.

Acknowledgements

We thank Angela Adamo and Søren Larsen for helpful comments and suggestions, David Cook for providing data for the dwarf galaxy sample and the anonymous referee for a constructive report which improved the paper. JP, NB and CU gratefully acknowledge funding from a European Research Council consolidator grant (ERC-CoG-646928-Multi-Pop). NB and RAC are Royal Society University Research Fellows. JMDK acknowledges funding from the German Research Foundation (DFG - Emmy Noether Research Group KR4801/1-1). JMDK and MRC acknowledge funding from the European Research Council (ERC-StG-714907, MUSTANG). MRC is supported by a PhD Fellowship from the International Max Planck Research School at the University of Heidelberg (IMPRS-HD). 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.

Appendix A Relationships between the CFE, pressure and $\Sigma_{\mathrm{SFR}}$

In Fig. 13 we compare the relationships between the CFE, natal gas pressure and $\Sigma_{\mathrm{SFR}}$ for the simulations. The top panel of the figure shows the CFE as a function of the natal gas pressure, which is the fundamental relationship underpinning the models. The galaxies closely trace the relation between CFE and pressure used at the particle level in the simulations (dash-dotted line in the figure, see Pfeffer et al., 2018), with a small amount of scatter caused by pressure fluctuations within the galaxies. In this parameter space, metallicity plays no role and galaxies of all metallicities follow the relation.

In the middle panel of Fig. 13 we again compare the CFE as a function of $\Sigma_{\mathrm{SFR}}$ as in Fig. 2, but with symbols coloured by metallicity. The scatter in this parameter space is significantly larger than for CFE versus pressure (top panel), demonstrating that scatter in the relation between $\Sigma_{\mathrm{SFR}}$ and pressure largely determines the scatter in the CFE-$\Sigma_{\mathrm{SFR}}$ relation. We show this explicitly in the bottom panel of Fig. 13, where we compare the natal gas pressure with $\Sigma_{\mathrm{SFR}}$ for the galaxies. The black dashed line in this panel shows the expected relationship, which describes the star formation relation used in the EAGLE model (Schaye & Dalla Vecchia, 2008; Schaye et al., 2015). At $\Sigma_{\mathrm{SFR}}\gtrsim 10^{-2}\,\mathrm{M}_{\sun}\,\mathrm{yr}^{-1}\,\mathrm{kpc}^{-2}$, most galaxies follow the Schaye & Dalla Vecchia (2008) relation. However at $\Sigma_{\mathrm{SFR}}\lesssim 10^{-2.5}\,\mathrm{M}_{\sun}\,\mathrm{yr}^{-1}\,\mathrm{kpc}^{-2}$, the galaxies are systematically offset from the relation (and this is also the case for a smaller fraction of galaxies at higher $\Sigma_{\mathrm{SFR}}$). This offset does not depend on the definition for $R_{\mathrm{lim}}$ (Section 2.2) and occurs for both limits of $1.5R_{1/2}$ and the radius containing 68 per cent of the recent star formation. This deviation is largely caused by highly concentrated or substructured star formation, such that star formation occurs in a much smaller area compared to the area for which $\Sigma_{\mathrm{SFR}}$ is calculated. The inverse is not the case (too high pressures for a given $\Sigma_{\mathrm{SFR}}$), since smaller apertures can always be chosen to bring the galaxies into better agreement with the expected $P$-$\Sigma_{\mathrm{SFR}}$ relation (at the expense of star particle and cluster numbers). Such an effect can be mitigated by calculating a mass-weighted surface density (see Johnson et al., 2016). In Fig. 14, we show the natal gas pressure compared with the mass-weighted $\Sigma_{\mathrm{SFR}}$ for the simulated galaxies (computed within $R_{\mathrm{lim}}$, as for Fig. 13). We calculate the mass-weighted $\Sigma_{\mathrm{SFR}}$ in a grid of $0.7\times 0.7\,\mathrm{kpc}^{2}$ regions (large enough such that incomplete sampling of star-forming regions is not important, Kruijssen & Longmore, 2014), with each region weighted by the mass of young ($<300\,\mathrm{Myr}$) stars. Using a mass-weighted $\Sigma_{\mathrm{SFR}}$ brings the low pressure ($P/k\lesssim 5\times 10^{3}\,{\rm K}\,{\rm cm}^{-3}$) galaxies in line with the expected $P$-$\Sigma_{\mathrm{SFR}}$ relation (with the exception of low-metallicity galaxies, see below), demonstrating the effectiveness of the method in accounting for spatially non-uniform star formation in the galaxies.

A secondary effect causing offset in the pressure-$\Sigma_{\mathrm{SFR}}$ relation is due to the metallicity of star formation. The EAGLE model includes a metallicity-dependent density threshold for star formation, such that star formation must occur at higher densities at lower metallicity (Schaye et al., 2015). This threshold is included to model the effect of the thermogravitational collapse of warm, photoionized interstellar gas into a cold, dense phase, which is expected to occur at lower densities and pressures in metal-rich gas (Schaye, 2004). This higher density threshold at lower metallicities results (through the lower density limit for star formation imposed by the polytropic equation-of-state implemented at high gas densities) in higher pressures of star formation at a given $\Sigma_{\mathrm{SFR}}$, which is evident at $\Sigma_{\mathrm{SFR}}\lesssim 5\times 10^{-3}\,\mathrm{M}_{\sun}\,\mathrm{yr}^{-1}\,\mathrm{kpc}^{-2}$ in the lower panel of Fig. 13 and Fig. 14. The threshold is also responsible for the apparent pressure floor in the figures (at $P/k\sim 10^{3}\,{\rm K}\,{\rm cm}^{-3}$). Note that this effect of increasing pressure with metallicity is not expected to occur at high $\Sigma_{\mathrm{SFR}}$, since the densities of star-forming gas in this regime are well above the metallicity-dependent threshold. The effect is also more apparent at higher redshifts in the simulations, where the metallicities of star-forming gas in the galaxies are lower.

In the bottom panel of Fig. 13, we also show the pressure-$\Sigma_{\mathrm{SFR}}$ adopted in the fiducial Kruijssen (2012) model (Krumholz & McKee, 2005). This relation is offset to lower $\Sigma_{\mathrm{SFR}}$ by $\approx 0.6$ dex when compared to the relation used in the EAGLE model. In the middle panel of Fig. 13, we also show the Kruijssen (2012) relation shifted to higher $\Sigma_{\mathrm{SFR}}$ to account for this offset (black dashed line). The fiducial (grey dashed line) and shifted (black dashed line) relations for the Kruijssen (2012) model thus give an indication of how the uncertainty in the pressure-$\Sigma_{\mathrm{SFR}}$ relation results in an uncertainty in the CFE-$\Sigma_{\mathrm{SFR}}$ relation.

Appendix B Truncations of initial and final cluster mass functions

In Fig. 15, we test how our fits of the cluster mass function truncation are affected by cluster mass loss in old, evolved cluster populations. We compare the upper exponential truncation mass ($M_{\mathrm{c,\ast}}$) resulting from Schechter (1976) function fits to the initial and evolved cluster populations in the simulated galaxies. For fits to the initial cluster masses, we fix the lower-mass power-law index to $\beta=-2$ (the initial value set in the simulations). For fits to the final cluster masses, we use a prior on the index $\beta$ between $-3$ and $-0.5$. The figure only shows galaxies where the most massive cluster is more massive than $M_{\mathrm{c,\ast}}$ (the median of the posterior distribution), for both initial and final cluster mass fits (see the discussion in Section 3.2). Overall, we find good agreement in $M_{\mathrm{c,\ast}}$ between fitting initial and evolved cluster masses, with potentially an offset to higher initial $M_{\mathrm{c,\ast}}$. However, a small offset to higher initial $M_{\mathrm{c,\ast}}$ is expected simply due to stellar-evolutionary mass loss (a factor $\sim 1.3$-$1.4$ for clusters with ages 100-300 Myr, or $0.1$ dex, based on the mass loss in the simulations). A number of galaxies have large errorbars for the fits to the final cluster masses due to the degeneracy between the truncation mass $M_{\mathrm{c,\ast}}$ and index $\beta$.

Appendix C Relation of $\phi_{P}$ to galactic properties

In Section 3.2 we find that, at a given $\Sigma_{\mathrm{SFR}}$, $M_{\mathrm{c,\ast}}$ is elevated above the $z=0$ relation at higher redshifts. We suggest that this effect may be caused by the higher gas fractions in galaxies at higher redshifts, through the parameter $\phi_{P}$ (Krumholz & McKee, 2005), and therefore in this section we directly test that suggestion.

In Fig. 16, we compare the median local $\phi_{P}$ for recently formed stars ($<300\,\mathrm{Myr}$) with $\Sigma_{\mathrm{SFR}}$ for all galaxies in Figs. 4 and 5. We find that $\phi_{P}$ does not directly correlate with $\Sigma_{\mathrm{SFR}}$, but instead with the gas fraction, which we show in Fig. 17. Therefore, at high redshift or low $\Sigma_{\mathrm{SFR}}$, we typically find $\phi_{P}\approx 1$, since the galaxies have low (stellar) mass with high gas fractions. The median $\phi_{P}$ increases from $\phi_{P}\approx 1$ at $z=2$ to $\phi_{P}\approx 2.5$ at $z=0$, implying a factor of four increase in the Toomre mass (since $M_{\mathrm{T}}\propto\phi_{P}^{-1.5}$, see equations 6 and 7 in Pfeffer et al., 2018), and thus in $M_{\mathrm{c,\ast}}$ (at the same pressure and epicyclic frequency $\kappa$). We show this effect in Fig. 5 as dashed lines, by scaling the fit to observed local galaxies (Johnson et al., 2016) by the median $\phi_{P}$ at each redshift. Therefore, the decrease of the typical $\phi_{P}$ with increasing redshift explains the elevation of $M_{\mathrm{c,\ast}}$ at higher redshifts found for the simulations.