Origins of scaling relations of globular cluster systems

TL;DR

This paper presents a detailed analysis of the physical origins and evolution of the scaling relations of globular cluster systems, using an analytic model that matches observed relations across a wide range of galaxy masses.

Contribution

It provides new insights into the evolution of the GC-to-halo mass ratio, the role of accreted satellites, and the formation history of metal-poor and metal-rich clusters.

Findings

The GC mass to halo mass ratio is nearly constant over redshift but evolves in normalization.

The shape of the stellar mass-halo mass relation influences the GC specific frequency.

Accreted satellites contribute increasingly to GC systems in more massive galaxies.

Abstract

Globular cluster (GC) systems demonstrate tight scaling relations with the properties of their host galaxies. In previous work, we developed an analytic model for GC formation in a cosmological context and showed that it matches nearly all of the observed scaling relations across 4 orders of magnitude in host galaxy mass. Motivated by the success of this model, we investigate in detail the physical origins and evolution of these scaling relations. The ratio of the combined mass in GCs $M_{\rm GC}$ to the host dark matter halo mass $M_h$ is nearly constant at all redshifts, but its normalization evolves by a factor of $\sim$10 from birth to $z=0$. The relation is steeper than linear at halo masses $M_h \lesssim 10^{11.5} M_{\odot}$, primarily due to non-linearity in the stellar mass-halo mass relation. The near constancy of the ratio $M_{\rm GC}/M_h$, combined with the shape of the stellar mass-halo mass relation, sets the characteristic $U-$shape of the GC specific frequency as a function of host galaxy mass. The contribution of accreted satellite galaxies to the buildup of GC systems is a strong function of the host galaxy mass, ranging from $\approx$0% at $M_h \approx 10^{11} M_{\odot}$ to 80% at $M_h \approx 10^{15} M_{\odot}$. The metal-poor clusters are significantly more likely to form ex-situ relative to the metal-rich clusters, but a substantial fraction of metal-poor clusters still form in-situ in lower mass galaxies. Similarly, the fraction of red clusters increases from $\approx 10$% at $M_h = 10^{11} M_{\odot}$ to $\approx 60$% at $M_h \approx 10^{13} M_{\odot}$, and flattens at higher $M_h$. Clusters formation occurs essentially continuously at high redshift, while at low redshift galactic mergers become increasingly important for cluster formation.

Figures

Tables

Table 1. Table 1: Accreted mass fractions for several bins of halo mass, shown in Fig. 6 . Column 1 gives the halo mass bin under consideration. Column 2 gives the interquartile range of the accreted mass fraction for all GCs. Columns 3 and 4 give the median accreted mass fractions for blue and red GCs, respectively. Columns 5 and 6 give the accreted stellar mass fractions for the field, from the Illustris simulation (Rodriguez-Gomez et al., 2016 ) and the UniverseMachine (Behroozi et al., 2018 ) .

Range of ${\log }_{10}{M}_h/M_{\odot }$	$f_{\mathrm{sat}}$ [IQR]	$f_{\mathrm{sat},\mathrm{blue}}$	$f_{\mathrm{sat},\mathrm{red}}$	$f_{\mathrm{sat},\mathrm{field}}$ [Illustris]	$f_{\mathrm{sat},\mathrm{field}}$ [UniverseMachine]
$11-11.5$	0.0 - 0.19	0.0	0.0	0.036	0.020
$11.5-12$	0.035 - 0.31	0.17	0.0	0.091	0.034
$12-12.5$	0.12 - 0.40	0.41	0.16	0.16	0.11
$12.5-13$	0.19 - 0.49	0.62	0.16	0.28	0.29
$13-13.5$	0.35 - 0.61	0.78	0.36	0.42	0.47
$13.5-14$	0.47 - 0.70	0.86	0.56	0.56	0.62
$14-14.5$	0.61 - 0.80	0.92	0.65	0.69	0.67

Equations

$$M_{\mathrm{GC}}\approx 3.4\times 10^{-5}{M}_{\mathrm{h}}.$$

$$R_{\mathrm{m}}\equiv \frac{M_{\mathrm{h},2}-M_{\mathrm{h},1}}{t_2-t_1}.\frac{1}{M_{\mathrm{h},1}},$$

$$M_{\mathrm{tot}}=1.8\times 10^{-4}{p}_2{M}_g,$$

$$\frac{M_g}{M_{\star }}(M_{\star },z)=0.35\times 3^{2.7}{\left(\frac{M_{\star }}{10^9{M}_{\odot }}\right)}^{-n_m(M_{\star })}{\left(\frac{1+z}{3}\right)}^{n_z(z)},$$

$$n_z$$

$$n_m$$

$$[Fe/H]={\log }_{10}\left[{\left(\frac{M_{\star }}{10^{10.5}{M}_{\odot }}\right)}^{0.35}(1+z{)}^{-0.9}\right].$$

$$\frac{dN}{dM}=M_0{M}^{-2}{e}^{-M/M_c},$$

$$\frac{dM}{dt}=-\frac{M}{t_{\mathrm{tid}}},$$

$$t_{\mathrm{tid}}(M)=5\mathrm{Gyr}{\left(\frac{M(t)}{2\times 10^5{M}_{\odot }}\right)}^{2/3}\left(\frac{P}{0.5}\right),$$

$$M(t)=M(t=0){\left[1-\frac{2}{3}\frac{t}{t_{\mathrm{tid}}(t=0)}\right]}^{3/2}\times f_{\mathrm{lost},\mathrm{se}}(t),$$

$$T_N\propto \frac{M_{\mathrm{GC}}}{M_{\mathrm{h}}}\frac{M_{\mathrm{h}}}{M_{\star }}.$$

$$\mathcal{M}\equiv \frac{1}{N_h}\sum _h{\left(\frac{M_{\mathrm{GC}}(z=0)}{M_{\mathrm{GC},\mathrm{obs}}(M_{\mathrm{h}})}-1\right)}^2+\frac{1}{N_h}\sum _h{\left(\frac{0.58}{{\sigma }_{Z,h}}\right)}^2+\frac{1}{G_M}+\frac{2}{G_Z},$$

Full text

††pubyear: 2019††pagerange: Origins of scaling relations of globular cluster systems–Origins of scaling relations of globular cluster systems

1 Introduction

The combination of their high masses, small sizes, and old ages means that globular clusters (GCs) in the local universe can be used as fossils of the extreme, high-redshift environments in which they formed. Although the upcoming James Webb Space Telescope and future thirty-meter class telescopes are poised to bring direct observations of GC formation at high-redshift (Renzini, 2017; Boylan-Kolchin, 2018; Zick et al., 2018; Vanzella et al., 2017, 2019), current observations of GCs are limited to within 200 Mpc (Harris et al., 2014). Inside this range, however, there is now a wealth of information on GC systems (e.g., Huchra et al., 1991; Harris, 1996; Beasley & Sharples, 2000; Beasley et al., 2000; Peng et al., 2006; Leaman et al., 2013; Brodie et al., 2014; Harris et al., 2014; Harris et al., 2016; Harris et al., 2017a; Forbes et al., 2018a). Despite the diversity of systems studied, observations have revealed that the ensemble properties of GC systems consistently correlate with the properties of their host galaxies. Below, we briefly summarize the main properties of GC systems relevant to this work.

The combined mass in GCs $M_{\mathrm{GC}}$ has been shown to scale almost linearly with the total host galaxy mass $M_{\mathrm{h}}$ (including the dark matter halo) over at least 4 orders of magnitude, with a normalization (Spitler & Forbes, 2009; Hudson et al., 2014; Harris et al., 2017b; Forbes et al., 2017; Forbes et al., 2018b):

[TABLE]

On the other hand, the specific frequency of GCs (the number of GCs divided by the galaxy stellar mass or luminosity) is highly non-linear and shows a characteristic $U$-shape, which is minimized for galaxies of stellar mass $M_{\star}\sim 10^{10}\,M_{\odot}$. Moreover, within an individual galaxy, the specific frequency also decreases with metallicity (Harris & Harris, 2002; Beasley et al., 2008; Lamers et al., 2017).

In contrast to the relatively universal GC mass function, the metallicity distribution function (MDF) of GC systems varies greatly from galaxy to galaxy. Because of the difficulty in acquiring high-resolution spectroscopy of extragalactic GCs, these MDFs are typically derived from integrated photometry combined with an empirically calibrated colour-metallicity relation (Peng et al., 2006; Harris et al., 2006; Usher et al., 2012). The mean and dispersion of GC system MDFs have also been found to scale with the mass of the host galaxy over a large range in galaxy mass (Peng et al., 2006). Qualitatively, many galaxies have clearly multi-modal MDF (e.g., Gebhardt & Kissler-Patig, 1999; Peng et al., 2006; Brodie & Strader, 2006; Brodie et al., 2012), but the most massive galaxies typically show broad, single-peaked distributions (Harris et al., 2014; Harris et al., 2016; Harris et al., 2017a).

The discovery of multi-modality in GC system MDF has historically led to the division of GCs into the so-called metal-poor “blue” and metal-rich “red” subpopulations, with a typical dividing line at $\mathrm{[Fe/H]}\approx-1$. These two subpopulations also differ in their spatial and kinematic properties. The blue GCs are typically on radially extended, eccentric orbits supported by random motions, while the red GCs follow more closely the field star distribution and kinematics (Zinn, 1985; Strader et al., 2011; Forbes et al., 2011; Pota et al., 2013; Durrell et al., 2014; Schuberth et al., 2010). Motivated by these observations, many authors have since suggested that the blue GCs formed ex-situ, in now-disrupted satellite galaxies, while the red GCs formed primarily in-situ (e.g., Cote et al., 1998; Katz & Ricotti, 2014). On the theoretical side, recent efforts have sought to build unified models of GC formation as a natural outcome of normal star formation occurring at the highest densities and pressures Kravtsov & Gnedin (2005); Muratov & Gnedin (2010); Tonini (2013); Li & Gnedin (2014); Kruijssen (2015); Li et al. (2017); Pfeffer et al. (2018); Choksi et al. (2018), rather than invoking distinct formation physics for each of the two subpopulations (as proposed by, e.g., Forbes et al., 1997; Cote et al., 1998; Beasley et al., 2002; Strader et al., 2005; Griffen et al., 2010).

Observations of ubiquitous young massive star clusters by the Hubble Space Telescope in nearby interacting and starbursting galaxies, such as the spectacular Antennae galaxies (e.g., Whitmore et al., 1993), have solidified the idea that mergers of gas-rich galaxies may trigger formation of star clusters destined to evolve into globular clusters. Galaxy mergers can pressurize the ISM and provide tidal torques that ultimately drive large-scale inflows of cold gas (Ashman & Zepf, 1992; Kravtsov & Gnedin, 2005; Bournaud et al., 2008; Renaud et al., 2015; Li et al., 2017), increasing the bound cluster formation efficiency (Kruijssen, 2012).

In Choksi et al. (2018), we presented an analytic model for the formation and evolution of GC systems along these lines. In this model, periods of rapid accretion onto dark matter halos are assumed to trigger cluster formation. The properties of the cluster population that form are set by the properties of the host galaxy, which are in turn set by empirical scaling relations derived from the observational literature. We found that GCs form over a wide range of redshifts and host galaxy masses, but the majority of cluster formation still occurs before the peak of field star formation. Furthermore, we demonstrated that this simple model simultaneously matches nearly all of the scaling relations obeyed by GC systems.

In Choksi & Gnedin (2019) we improved on the physical realism of this model by updating the adopted cluster initial mass function from a pure power-law to a Schechter (1976) function $dN/dM\propto M^{-2}e^{-M/M_{c}}$ which exponentially truncates the formation of the most massive clusters. This change was motivated by numerous observations of the initial mass function of young clusters in nearby galaxies (e.g., Gieles et al., 2006; Larsen, 2009; Bastian, 2008; Adamo et al., 2015; Johnson et al., 2017) and by galaxy formation simulations that model cluster formation (Li et al., 2017; Li et al., 2018). After experimenting with various values of the cutoff mass, we determined that $M_{c}\approx 10^{6.5}-10^{7}\,M_{\odot}$ produced results that were most consistent with present-day GC scaling relations. Lower values of $M_{c}$ prevent the formation of very massive GCs found in massive elliptical galaxies (Harris et al., 2014), while higher values are disfavoured by modeling the high-mass end of the present-day GC mass function (Johnson et al., 2017). Section 2 provides more detail on the model setup.

In this work, we investigate the origin of the scaling relations presented in Choksi et al. (2018). We begin in Section 3 by examining the evolution over time of the GC system mass - halo mass ($M_{\mathrm{GC}}-M_{\mathrm{h}}$) relation. We also discuss the origins of non-linearity in the relation for dwarf galaxies and the implications for the specific frequencies of GCs. In Section 4 we investigate the contribution of accreted satellite galaxies and their GC systems to the properties of present-day GC systems. In Section 5 we discuss the role of galaxy mergers in GC formation. Finally, in Section 6 we discuss the implications of our results and summarize our conclusions.

2 Methodology

Below, we briefly summarize the main components of our cluster formation model. For further details and justification regarding the choices of equations and parameters, we refer readers to Choksi et al. (2018) and Choksi & Gnedin (2019).

2.1 Summary of cluster formation model

We trigger cluster formation when the specific mass accretion rate $R_{\rm m}$ onto a dark matter halo exceeds a threshold value $p_{3}$ between consecutive outputs of our adopted collisionless cosmological simulation. Specifically, for a halo of mass $M_{\rm h,2}$ at time $t_{2}$ and its progenitor of mass $M_{\rm h,1}$ at time $t_{1}$, we calculate the specific mass accretion rate as:

[TABLE]

and trigger cluster formation for $R_{\rm m}>p_{3}$. When cluster formation is triggered, we form a population of clusters of combined mass $M_{\mathrm{tot}}$:

[TABLE]

where $M_{g}$ is the cold gas mass in the galaxy and the normalization of equation (missing) 2 is motivated by the cosmological hydrodynamic simulations of Kravtsov & Gnedin (2005). The values of $p_{2}$ and $p_{3}$ are taken to be free parameters111The current form of our model has only two adjustable parameters, but we keep the labels $p_{2}$ and $p_{3}$ for consistency with previous versions of the model., whose values are constrained using a wide array of observational data on GC system metallicities and masses (Harris, 1996; Côté et al., 2006; Peng et al., 2006; Harris et al., 2014; Harris et al., 2016; Harris et al., 2017a).

The cold gas fraction is parameterized as a function of the stellar mass $M_{\star}$ and redshift $z$ as:

[TABLE]

where $n_{z}$ and $n_{m}$ are given by:

[TABLE]

The host galaxy stellar mass is increased self-consistently using a modified version of the semi-empirical stellar mass-halo mass relation derived from forward modeling by Behroozi et al. (2013), and extrapolated at high-redshift ($z\gtrsim 8$). The stellar mass also sets the galaxy metallicity via an empirical galaxy mass-metallicity relation:

[TABLE]

Individual cluster metallicities are set by the metallicity of the host galaxy in which they form, with an additional scatter of 0.3 dex. Clusters are drawn from a cluster initial mass function of the form:

[TABLE]

where $M_{0}$ is an overall normalization factor and $M_{c}$ is the characteristic truncation mass. Our method for sampling the cluster initial mass function is the “optimal sampling” method of Schulz et al. (2015).

In Choksi & Gnedin (2019) we analyzed the effects of different values of $M_{c}$ on GC system scaling relations. However, because it is not the focus of this work, here we present results from only one model of fixed $M_{c}=10^{7}\,M_{\odot}$. This value of $M_{c}$ is consistent with theoretical expectations and inferences from observations of the local GC mass functions, and also robustly reproduces the various GC system scaling relations (Reina-Campos & Kruijssen, 2017; Jordán et al., 2007; Johnson et al., 2017; Choksi & Gnedin, 2019). All results presented in this work are qualitatively robust to variations in $M_{c}$. For the $M_{c}=10^{7}\,M_{\odot}$ model we adopt in this work, the best-fit values of the free parameters are $p_{2}=8.8$, $p_{3}=0.58\rm\,Gyr^{-1}$.

GCs lose mass gradually due to two-body relaxation and tidal stripping. Because our model contains no spatial information on cluster orbits, we apply a spatially-averaged dynamical disruption prescription following Gnedin et al. (2014):

[TABLE]

where the disruption timescale $t_{\mathrm{tid}}$ was calibrated using direct $N-$body simulations (Gieles & Baumgardt, 2008):

[TABLE]

where $P$ is a normalized period of rotation around the galactic center. In Choksi & Gnedin (2019) we showed that adopting a constant value of $P=0.5$ reproduces well the present-day GC mass function. Integrating equation (missing) 6 yields the mass as a function of time $t$ since cluster birth:

[TABLE]

where we apply a fractional mass loss due to stellar evolution $f_{\rm lost,se}(t)$ as calculated by Prieto & Gnedin (2008).

3 The globular cluster system mass - halo mass relation

In this section we investigate the origin and evolution of the relationship between the combined mass in GCs and the host halo mass. Throughout this section and next, we distinguish between the central galaxy and its satellites using the “main progenitor branch” (MPB) tag of the adopted dark matter merger trees. The MPB is defined as the branch along the merger tree with the largest integrated mass history (for details, see De Lucia & Blaizot, 2007). It can be thought of as a “current” halo mass at each redshift. Merger trees are taken from the collisionless run of the Illustris simulation and constructed using the sublink algorithm (Springel et al., 2001; Vogelsberger et al., 2014).

Fig. 1 shows the $M_{\mathrm{GC}}-M_{\mathrm{h}}$ relation at three distinct epochs, $z=10$, 3, and 0. We include in the calculation of $M_{\mathrm{GC}}$ only the clusters that have already been accreted onto the MPB, so that $M_{\mathrm{GC}}$ represents the observable mass in clusters in the galaxy at a given epoch, including the effect of disruption and stellar evolution. $M_{\mathrm{h}}$ is the mass of the MPB halo of each merger tree. We find that the nearly linear relationship between $M_{\mathrm{GC}}$ and $M_{\mathrm{h}}$ holds across redshifts, but with a normalization $M_{\mathrm{GC}}/M_{\mathrm{h}}$ that decreases monotonically with cosmic time.

When cluster formation is triggered, we form a population of clusters of combined mass $M_{\mathrm{tot}}\sim 1.6\times 10^{-3}M_{g}$. At $z\sim 10$, near the beginning of cluster formation, most galaxies are very gas-rich, and $M_{g}\sim f_{b}M_{\mathrm{h}}$. Each cluster formation event then produces a mass in clusters $M_{\mathrm{tot}}\sim 3\times 10^{-4}M_{\mathrm{h}}$. There are several outputs of the Illustris simulation at $z>10$, which may combine the young cluster populations to reach $M_{\mathrm{GC}}\sim 10^{-3}M_{\mathrm{h}}$. This is the upper envelope of the normalization seen at high redshifts in Fig. 1.

At high redshift, continuous cluster formation and accretion of external GC systems during galaxy mergers, combined with concurrent growth of the host halo, causes points to move diagonally in the $M_{\mathrm{GC}}-M_{\mathrm{h}}$ plane. The normalization decreases only by a factor of $\sim 2$ from $z=10$ to $z=3$.

On the other hand, we find that between $z=3$ and $z=0$ the normalization decreases by a much larger factor of $\sim 10$: from $M_{\mathrm{GC}}\approx 3\times 10^{-4}M_{\mathrm{h}}$ at $z=3$ to $M_{\mathrm{GC}}\approx 3\times 10^{-5}M_{\mathrm{h}}$ at $z=0$. Over this redshift range, the GC system mass can decrease because of dynamical disruption, or increase because of late cluster formation, but on average it remains nearly constant. In contrast, the halo mass only increases with time and this dominates the evolution of the normalization at $z\lesssim 3$ and causes the points in Fig. 1 to move rightward in the $M_{\mathrm{GC}}-M_{\mathrm{h}}$ plane. Furthermore, because the logarithmic accretion rates onto dark matter halos depend only weakly on halo mass, $\dot{M_{\mathrm{h}}}/M_{\mathrm{h}}\propto M_{\mathrm{h}}^{0.13}$ (McBride et al., 2009), the points shift rightward about the same amount and therefore the shape of the $M_{\mathrm{GC}}-M_{\mathrm{h}}$ relation is preserved until $z=0$.

However, our model also predicts deviation from a purely linear relation. At halo masses below $M_{\mathrm{h}}\sim 10^{11.5}\,M_{\odot}$ the predicted $M_{\mathrm{GC}}-M_{\mathrm{h}}$ scaling falls off faster than linearly. Moreover, this behaviour is already present at $z\sim 10$, meaning that it is imprinted at birth and is not due to how the GC system evolves. The shape of the relation between the cold gas mass and halo mass therefore is the dominant factor in setting the model’s non-linearity. The $M_{g}-M_{\mathrm{h}}$ relation is set indirectly, via the $M_{g}-M_{\star}$ relation in equation (missing) 3 and the $M_{\star}-M_{\mathrm{h}}$ relation of Behroozi et al. (2013). The latter of these is highly nonlinear and dominates the nonlinear behaviour seen in the model. Therefore, the observed non-linearity in Fig. 1 can be used to place a constraint on the shape of the stellar mass-halo mass relation at high redshift.

To quantify the degree of non-linearity, we perform a linear fit to the model $M_{\mathrm{GC}}-M_{\mathrm{h}}$ relation at each redshift and calculate the median deviation (in perpendicular log-mass distance) of the model points from the best-fit line. At $z>3$, this deviation remains at small values 0.13-0.14 dex, barely distinguishable from the linear relation. By $z=0$ the median deviation increases to 0.19 dex, and the relation can be considered mildly non-linear.

At the largest galaxy masses, the $M_{\star}-M_{\mathrm{h}}$ relation (Behroozi et al., 2013) is significantly more non-linear than the $M_{\mathrm{GC}}-M_{\mathrm{h}}$ relation. Giant early-type galaxies are believed to assemble by consuming a large number of satellite galaxies, thus combining both their field stellar populations and GC systems. Why does this process result in visibly different scaling relations for $M_{\star}$ and $M_{\mathrm{GC}}$? We return to this point in next section where we investigate the contribution of satellite galaxies to the overall GC system.

Finally, although the non-linearity is already present at birth, we note that the lower average cluster masses in low-mass galaxies (Choksi et al., 2018), combined with the fact that low-mass clusters have larger fractional mass loss by equation (missing) 7, may exacerbate the steeper than linear fall-off in dwarf galaxies. We plan to investigate the GC systems of dwarf galaxies in more detail in future work.

A commonly used statistic in studying GC systems is the specific frequency $T_{N}\equiv N_{\rm GC}/(M_{\star}/10^{9}\,M_{\odot})$. This specific frequency is proportional to

[TABLE]

The first ratio in the above expression, $M_{\mathrm{GC}}/M_{\mathrm{h}}$, is nearly constant. The ratio $M_{\mathrm{h}}/M_{\star}$, however, is a strong function of $M_{\star}$ and is minimized for $M_{\star}\approx 10^{10.2}\,M_{\odot}$ (e.g., Behroozi et al., 2013). Thus, in our model the $U-$shape of the specific frequency as a function of galaxy mass is predominantly set by the shape of the stellar-mass halo mass relation (see also El-Badry et al. 2019). Fig. 2 shows that our model prediction for $T_{N}$ is consistent with data from the Virgo Cluster Survey (as presented in Peng et al. 2008). This is expected given that the model matches the $M_{\mathrm{GC}}-M_{\mathrm{h}}$ relation. At low halo masses, the ratio $M_{\mathrm{GC}}/M_{\mathrm{h}}$ is no longer constant and its variation with host mass will begin to affect the shape of the $T_{N}$ relation.

4 Effects of satellite systems

In this section we examine the impact of accreted satellite systems on the buildup of present-day GC systems. As discussed in Section 3, we identify the central galaxy and its satellites at each redshift using the “main progenitor branch” tag of the adopted dark matter merger trees.

Fig. 3 illustrates the effect of satellite GC systems on the $M_{\mathrm{GC}}-M_{\mathrm{h}}$ relation. It shows the fiducial model prediction with and without contribution of GCs formed in satellite systems. At low halo masses, the two cases are indistinguishable. However, the difference becomes noticeable in the most massive (group or cluster sized) halos with $M_{\mathrm{h}}\gtrsim 10^{13}\,M_{\odot}$. Without the contribution of accreted GCs, the GC system mass is too low by up to 0.3 dex.

Fig. 4 demonstrates the more drastic impact of satellite accretion on the MDF of GC systems. Without the contribution of accreted systems, the mean metallicity would continue to rise with halo mass, rather than plateauing at $\mathrm{[Fe/H]}\approx-0.8$. The typical mean metallicity of GC systems in the most massive halos with $M_{\mathrm{h}}\approx 10^{14.5}\,M_{\odot}$ would increase by about 0.5 dex relative to the observations and fiducial model.

A related quantity is the dispersion in metallicities of a GC system, $\sigma_{Z}$ (the width of the MDF) plotted in Fig. 5. In Choksi et al. (2018) we showed that the weak scaling of $\sigma_{Z}$ with $M_{\mathrm{h}}$ is a robust prediction of our model. However, when we exclude accreted GCs, we find a break in this relation at $M_{\mathrm{h}}\sim 10^{12.5}$, such that $\sigma_{Z}$ begins to decrease with halo mass. In the most massive galaxies, the width of the metallicity distribution is far too narrow without the contribution of accreted satellite systems. In lower mass halos below $M_{\mathrm{h}}\approx 10^{12.5}\,M_{\odot}$, $\sigma_{Z}$ still scales weakly with the halo mass because the main branch dominates cluster formation events.

The results visualized in Figs. 3-5 can be understood by analyzing Fig. 6. The shaded grey region in that figure shows the interquartile range for the accreted GC mass fraction for model halos. In low mass halos, in-situ GC formation dominates, while in giant halos the majority of GCs are formed ex-situ. This explains the increasing discrepancy between the “All” and “MPB Only” curves at high halo masses in the previous figures. Because satellites (and by extension, their GCs, according to the adopted galaxy mass-metallicity relation) will tend to have lower metallicities, excluding their contribution leads to a higher metallicity.

We note that the fraction of clusters by number formed in satellites follows a very similar trend for red clusters, but is always somewhat higher for blue clusters. This is because the mean mass of blue clusters formed in satellites is lower than that formed in the central galaxy by about a factor of 2. This is due to the fact that satellites have smaller gas reservoirs, which prevents sampling of the high-mass end of the cluster initial mass function (for a more detailed discussion of this effect see Choksi & Gnedin 2019, Appendix A). For a Milky Way mass halo with $M_{\mathrm{h}}\approx 10^{12}\,M_{\odot}$, we expect only $\approx$30% of GCs to have formed ex-situ.

The figure also shows two estimates of the accreted stellar mass fraction for all (“field”) stars. The dash-dot brown curve shows the semi-empirical result derived via forward modeling from the UniverseMachine model of Behroozi et al. (2018), while the dashed brown curve shows the prediction from the Illustris cosmological simulation as analyzed by Rodriguez-Gomez et al. (2016). These two curves can be compared to the shaded black region, showing the analogous result for GCs. All three curves show a similar normalization and slope for very massive galaxies, demonstrating that GCs indeed trace well the overall buildup of field stars in galaxies. However, in Milky Way-sized galaxies there remains a a factor of 2-3 difference in the median trend. Numerical values of the plotted fractions are given for reference in Table 1.

Fig. 6 also shows the median accreted GC mass fraction separately for the red and blue populations: $M_{\rm GC,blue,sat}/M_{\rm GC,blue}$ and $M_{\rm GC,red,sat}/M_{\rm GC,red}$. Clusters are identified by applying the Gaussian Mixture Modeling (GMM) algorithm of Muratov & Gnedin (2010) to the MDF of each GC system. We take the metallicity where the magnitude of two Gaussians are equal as the dividing line between blue and red clusters. This metallicity is typically located in the range $\mathrm{[Fe/H]}\approx-0.7$ to $-1.2$. Therefore, in systems with either too few clusters to reliably apply GMM or strongly unimodal MDFs, we continue to use a fixed cutoff $\mathrm{[Fe/H]}=-1.0$ to differentiate red and blue clusters.

As expected, we find that the blue clusters are systematically more likely to form ex-situ than red clusters. However, for a Milky Way mass galaxy we still predict a majority – roughly 70% – of the mass in blue clusters to have formed in-situ. On the other hand, in the most massive galaxies we predict the majority of mass of red clusters to have also formed ex-situ, due to the strong dependence of the accreted mass fraction with host halo mass.

To better understand the relative importance of the red and blue populations as a function of galaxy mass, we show in Fig. 7 the fractions (by number): $f_{\mathrm{red}}\equiv N_{\rm GC,red}/N_{\rm GC}$ and $f_{\rm blue}\equiv N_{\rm GC,blue}/N_{\rm GC}$. The model predicts $f_{\mathrm{red}}$ to increase with halo mass. In dwarf galaxies with $M_{\mathrm{h}}\approx 10^{11}\,M_{\odot}$, red clusters are highly subdominant, with $f_{\mathrm{red}}\approx 10\%$. On the other hand, in group and cluster environments with $M_{\mathrm{h}}\gtrsim 10^{13}\,M_{\odot}$, the model predicts more red than blue clusters, with $f_{\mathrm{red}}$ reaching an asymptotic value of about 60%. We note that El-Badry et al. (2019) also predicted the fraction of red clusters as a function of halo mass (see their Fig. 9). It is reassuring that despite the many differences in their inputs, both models predict a very similar behaviour of $f_{\mathrm{red}}$ with halo mass, including an asymptote at $f_{\mathrm{red}}\approx 60\%$ for $M_{\mathrm{h}}\gtrsim 10^{13}\,M_{\odot}$.

We have also overplotted in Fig. 7 upper limits (triangles) from the Virgo Cluster Survey (VCS; Peng et al., 2006) over the mass range $M_{\mathrm{h}}\approx 10^{11}-10^{13}\,M_{\odot}$. These points represent upper limits because HST covered only the central parts of their host galaxies, where red clusters preferentially reside. At $M_{\mathrm{h}}\gtrsim 10^{13}\,M_{\odot}$, we also show data from the recent HST survey of brightest cluster galaxies (Harris et al., 2014). Finally, we include data from the Milky Way and M31 (Harris, 1996; Huxor et al., 2014). The data are qualitatively consistent with the predictions of our model. For a Milky Way mass halo in particular, our model predicts $f_{\mathrm{red}}\approx 30\%$, in excellent agreement with the observed fraction in the Galaxy (Harris, 1996).

5 The role of galaxy mergers in triggering GC formation

In this section we analyze the role of mergers in our model. We emphasize that our criterion for forming clusters is only that the specific accretion rate onto a dark matter halo $R_{m}\equiv\dot{M_{\mathrm{h}}}/M_{\mathrm{h}}$ exceeds a threshold value $p_{3}$ which is calibrated based on observations; we do not explicitly require any link to actual mergers of two distinct halos. Therefore, we can test what fraction of model clusters form in events that satisfy the criterion $R_{m}>p_{3}$ and also represent major mergers.

To identify mergers, for each halo in the merger tree we search for progenitors at the previous simulation snapshot. If more than one progenitor is identified, we label the event as a merger. However, the typical spacing between outputs of our adopted simulation is only $\approx$0.1 Gyr, yet the time for dark matter halos to merge (the interval between being accreted and disappearing from the halo catalog) can range from 0.1 to 1 Gyr. Therefore, the merger may be resolved over the course of multiple simulation snapshots. To account for the extended duration of mergers, we consider all the snapshots within some time interval $\Delta t$ before the disappearance of the satellite from the halo catalog as “undergoing” a merger. The exact value of $\Delta t$ is somewhat arbitrary, as the range of merger timescales in the Illustris simulation is quite broad, peaking at $\sim 200$ Myr but having an extended tail up to 1 Gyr. For clarity we adopt a fixed duration $\Delta t=200$ Myr. A different choice of $\Delta t$ would affect the inferred merger probability (for example, $\Delta t=400$ Myr increases it by a factor up to 1.6) but would not qualitatively alter our conclusions.

For each merger we calculate the ratio of masses of the satellite and central halo, $q<1$. Because the merger may be resolved over multiple snapshots, the infalling satellite may have experienced significant tidal stripping before the last time it appears in the halo catalog. Therefore, the mass ratio calculated at the satellite’s last snapshot is not representative of the extent to which the galactic potential was disturbed by the merger. The final recorded satellite mass is also strongly resolution-dependent; for a complete discussion of these effects, see Rodriguez-Gomez et al. (2015). Instead, we calculate $q$ at the snapshot corresponding to the maximum mass of the satellite, which occurs approximately when the satellite crosses the virial radius of the central halo.222We verified our merger-identification process by calculating the cumulative number of major mergers along the main branch, in bins of $z=0$ halo mass, and compared this to the result in Fig. 7 of Fakhouri et al. (2010) for the Millennium-II simulation. We found generally good agreement in the redshift evolution and normalization of the number of mergers, with differences of $\lesssim 30\%$. The discrepancy is likely attributable to differences in the adopted collisionless simulation and details of the halo finding and merger identification process. Overall, this comparison gives us confidence in the number of major mergers we detect.

Applying this procedure, we find that $\approx$20% of all clusters with initial masses above $2\times 10^{5}\,M_{\odot}$ form during major merger events with $q>\nicefrac{{1}}{{4}}$, with only a weak scaling with the $z=0$ host halo mass.

To better understand this result, we calculated at each simulation snapshot the ratio of the number of halos that satisfy both conditions $R_{m}>p_{3}$ and $q>\nicefrac{{1}}{{4}}$ to the number of halos with $R_{m}>p_{3}$. In essence, this ratio represents the probability that a halo which is forming clusters is also undergoing a major merger, $P(q>\nicefrac{{1}}{{4}}\,|\,R_{m}>p_{3})$. Our model assumes that the number of clusters formed in any event is proportional only to the cold gas mass in the central galaxy, and does not include the actual effects of the merger (increase of ISM pressure, rapid change of gravitational potential, etc.). Therefore, this probability can be taken as a simple estimate of the fraction of clusters expected to form during major mergers, which relies only on $\Lambda$CDM statistics and not on the details of our GC formation model. The solid curves in Fig. 8 demonstrate that this probability ranges from $\approx 5\%$ for halos with $M_{\mathrm{h}}\sim 10^{11}\,M_{\odot}$ to $\approx 30\%$ for the highest-mass halos, with only a weak redshift evolution. On the other hand, the probability that a randomly selected halo is undergoing a major merger at a given redshift, $P(q>\nicefrac{{1}}{{4}})$, actually decreases steeply towards low redshift, as seen in the dashed curves in Fig. 8 which begin to diverge significantly from the solid curves at $z<2.5$.

Why does the probability of a halo undergoing a merger decrease with time while the probability of major merger-induced cluster formation remains constant? This difference indicates that the probability of a halo exceeding the threshold for forming clusters $P(R_{m}>p_{3})$ also systematically decreases with time. In terms of the probabilities it can be written as $P(q>\nicefrac{{1}}{{4}})=P(q>\nicefrac{{1}}{{4}}\,|\,R_{m}>p_{3})\,P(R_{m}>p_{3})$. In Fig. 9 we illustrate the sharp decline of the halo accretion rates $R_{m}$ with cosmic time, consistent with standard expectations from $\Lambda$CDM (e.g., McBride et al., 2009; Behroozi & Silk, 2015; Rodríguez-Puebla et al., 2016). At high redshift, $z>3$, typical halo accretion rates significantly exceed the threshold $p_{3}$ and cluster formation can occur almost continuously. Of course for some halos occasional drops in $R_{m}$ below $p_{3}$ can occur (as demonstrated by the dotted 5th percentile curve) but these are rare. On the other hand, at $z\lesssim 1.5$ the entire region of the $R_{m}$ interquartile range lies below the cluster formation threshold. However, a few formation events do still occur (as demonstrated by the dotted 95th percentile curve). Thus low-redshift cluster formation proceeds far more discretely.

6 Discussion

6.1 The origin of the GC-halo mass relation

El-Badry et al. (2019) showed that the tight $M_{\mathrm{GC}}-M_{\mathrm{h}}$ relation does not necessarily imply such a relationship at formation. The action of repeated galaxy mergers, which combine GC systems, naturally leads to a tight $M_{\mathrm{GC}}-M_{\mathrm{h}}$ relation by $z=0$ simply due to the central limit theorem. However, this does not preclude the existence of an $M_{\mathrm{GC}}-M_{\mathrm{h}}$ relation at formation; it only makes it unnecessary to explain local observations.

As we showed in Section 3, such a correlation is still expected at high redshift in our model, due to the adopted relation for the cluster formation rate $M_{\mathrm{tot}}\propto M_{g}$. This proportionality was initially motivated by the results of early cosmological simulations by Kravtsov & Gnedin (2005). Those simulations identified dense giant gas clouds in high-redshift galaxies and applied a sub-grid model for the formation of globular clusters. The sum of masses of all clusters formed in a given event was roughly proportional to both the halo mass and the baryon (gas plus stars) mass. The relation to gas mass is more directly causal because star clusters are observed to form from dense cold gas. For this reason, the original version of our model (Muratov & Gnedin, 2010) adopted the proportionality of $M_{\mathrm{tot}}$ to the mass of neutral gas in the host galaxy. In Choksi et al. (2018) we revised it to the mass of molecular gas, which should be even more closely related to the formation rate of massive clusters. While observations indicate that the surface density of molecular gas is most directly proportional to the star formation rate, our model cannot incorporate this relation because it does not calculate the spatial structure of galaxies. We believe our three-step scaling procedure from $M_{\mathrm{h}}$ to $M_{\star}$ to $M_{g}$ to $M_{\mathrm{tot}}$ is an acceptable parameterization of the cluster formation rate.

Kruijssen (2015) predicted the $M_{\mathrm{GC}}-M_{\mathrm{h}}$ relation at $z=3$ to be nearly identical to the $z=0$ relation, in contrast to our prediction of a factor of $\sim$10 evolution (see Fig. 1). His model uses the cluster formation efficiency of Kruijssen (2012) and assumes all GC formation happens in a single epoch. Furthermore, in contrast to our explanation of the shape of the $T_{N}-M_{\star}$ relation, Kruijssen (2015) uses the shape of $T_{N}(M_{\star})$ to explain the $M_{\mathrm{GC}}-M_{\mathrm{h}}$ relation and derive $M_{\mathrm{GC}}/M_{\mathrm{h}}\propto M_{\star}^{0.1}$. However, in our model, it is the $M_{\mathrm{GC}}-M_{\mathrm{h}}$ relation that is more fundamental, because the mass that forms in clusters (indirectly) scales almost linearly with the halo mass.

Similar to our results, the fiducial model of El-Badry et al. (2019) predicts a steeper-than-linear fall-off of the $M_{\mathrm{GC}}-M_{\mathrm{h}}$ relation at low galaxy mass. In their model, it is due to the cosmic UV background shutting off gas accretion onto dwarf galaxies. This effect is incorporated implicitly in our model in the shapes of the $M_{\star}-M_{\mathrm{h}}$ and $M_{\star}-M_{g}$ relations. In addition, our model includes the effect of systematically lower cluster masses in dwarf galaxies (Choksi et al., 2018), which further steepens the relation. However, the scatter of the relation increases dramatically at the dwarf galaxy scale, which makes comparison with observations more challenging and requires a more detailed investigation.

6.2 Accretion of satellite GC systems

Most massive early-type galaxies in the local universe are believed to form in two phases: a short high-redshift phase of intense, concentrated in-situ star formation followed by late time accretion of satellite galaxies (e.g., Oser et al., 2010; Forbes et al., 2011; van Dokkum et al., 2015). Recently, Beasley et al. (2018) studied the GC system of the nearby “red nugget” galaxy NGC 1277, so called because it is believed to be a low-redshift analog of early massive galaxies, i.e., it has not yet undergone the phase of late-time satellite accretion. They find NGC 1277 to be nearly devoid of blue GCs ($f_{\rm blue}\approx 17\%$) and present it as a puzzle for theory. This observation can be explained by our model, which predicts for typical galaxies of NGC 1277 mass ($M_{\star}\approx 10^{11}\,M_{\odot}$) that $\approx 40\%$ of all clusters are blue and $\approx 70\%$ of these blue clusters are accreted (Figs. 6-7). Thus the blue fraction could drop to as low as 12% if a galaxy happened to miss accretion of all satellites. This is plausible for NGC 1277 because it may have fallen early into the Perseus cluster and been tidally limited. In fact, the observed fraction $f_{\rm blue}$ indicates that some satellites did contribute to this GC system.

The predicted range of the fraction of blue clusters may thus be used to constrain the history of galaxy assembly. For the same stellar mass, galaxies in group or cluster environments where tidal limitation prevents accretion of satellite galaxies, should have a lower blue fraction than those in the field.

6.3 Impact of galaxy mergers on the cluster formation rate

Li et al. (2017); Li et al. (2018) used cosmological simulations of a single Milky Way mass halo to study in detail the formation process of massive star clusters. While our model ultimately produces similar results for the average GC formation epochs and demographics as in these simulations, the role of major mergers differs. Li & Gnedin (2019) find about 75% of clusters with initial masses $M\geqslant 2\times 10^{5}\,M_{\odot}$ form within 200 Myr of three major mergers, about a factor of 4 higher than our model prediction. The difference between these two results may stem from several sources. Galaxy formation simulations include a variety of detailed physics relevant for cluster formation that is not captured in our model. In particular, they find that the cluster initial mass function extends to higher masses during major mergers, increasing the probability of forming massive GCs that would survive to the present day. In contrast, our model adopts a merger-independent cluster formation rate. Finally, there is significant scatter in galaxy assembly histories, even at fixed galaxy mass, and therefore a single simulated realization may differ from the ensemble average of our model.

6.4 Which model parameters disfavour mergers?

Our current best-fit range on the free parameters $p_{2}$ and $p_{3}$ simultaneously matches properties of both the masses and metallicities of GC systems. Increasing the threshold accretion rate for forming clusters ($p_{3}$) and the cluster formation rate parameter ($p_{2}$) would lead to a higher fraction of all cluster formation occurring during major mergers. Such values could still match the observed $M_{\mathrm{GC}}-M_{\mathrm{h}}$ relation, by having fewer cluster formation events but more mass formed in each event. However, we found that such a parameter combination is disfavoured in our model. We have performed parameter optimization by minimizing the merit function (Choksi et al., 2018):

[TABLE]

where the first term is the reduced $\chi^{2}$ of the $M_{\mathrm{GC}}-M_{\mathrm{h}}$ relation, the second term compares the observed average $\sigma_{Z}$ of 0.58 dex to the $\sigma_{Z}$ of model galaxies, and $G_{M}$ and $G_{Z}$ are the mass and metallicity “goodness” statistics, defined as the fraction of model-observed galaxy pairs that have GC system metallicity and mass distributions with an acceptable Kolmogorov-Smirnov statistic.

To understand in detail what aspect of the model determines this conclusion, we investigated the relative contributions of each term in the merit function over a wide range $p_{2}<20$ and $p_{3}<3\,$Gyr*-1*. We found that across almost the entire parameter space the variation of the $\chi^{2}$ and $G_{Z}$ terms dominates the gradient of the merit function. The $G_{M}$ and $\sigma_{Z}$ terms vary negligibly and could be eliminated for simplicity without loss of model accuracy.

High values of $p_{2}$ and $p_{3}$, which would lead to a higher fraction of cluster formation during major mergers, still result in a higher $\chi^{2}$ term than the fiducial model. At the same time, high values of $p_{2}$ are strongly disfavored by the $G_{Z}$ term: a higher normalization of the cluster formation rate would make the GC system MDF closer to the field star MDF, which is peaked at a significantly higher metallicity than the observed GCs.

Finally, we note that Fig. 9 also demonstrates that the model results are not very sensitive to the exact value of $p_{3}$. While shifting $p_{3}$ slightly may change results in detail, the median $R_{m}$ evolves by two orders of magnitude over cosmic time. Therefore, the main function of $p_{3}$ is instead to select out the general range of accretion rates $R_{m}$ at low redshift that are allowed to form clusters.

7 Summary

In this work, we explored the origins and buildup of many of the GC system-host galaxy scaling relations. Our main results are:

The ratio $M_{\mathrm{GC}}/M_{\mathrm{h}}\sim 10^{-4}-10^{-3}$ at birth, and evolves only mildly at $z>3$ due to continuous cluster formation. From $z\approx 3$ to $z=0$, the ratio decreases by roughly a factor of ten because the GC formation rate drops sharply while the halo mass continues to increase (Fig. 1). 2. 2.

The shape of the relation between GC specific frequency $N_{\rm GC}/M_{\star}$ and host galaxy stellar mass exhibits a characteristic $U-$shape, which is set by the near constancy of the ratio $M_{\mathrm{GC}}/M_{\mathrm{h}}$ and the peaked shape of the stellar mass-halo mass relation (Fig. 2). 3. 3.

The fraction of GC mass accreted from now-disrupted satellite galaxies increases monotonically with $M_{\mathrm{h}}$, ranging from a few percent at $M_{\mathrm{h}}\approx 10^{11}\,M_{\odot}$ to 80% at $M_{\mathrm{h}}\approx 10^{14.5}\,M_{\odot}$. These values are similar to the accreted fraction of “field” stars in giant galaxies, but exceed them in Milky Way-sized and dwarf galaxies (Fig. 6 and Table 1). Blue GCs are systematically more likely to form in satellite galaxies than red GCs. 4. 4.

Consequently, without the contribution of accreted GCs the mean metallicities of GC systems would be up to 0.5 dex too high and the metallicity dispersions 0.4 dex too low (Figs. 4-5). The combined mass in GCs would also be up to 0.3 dex too low (Fig. 3). 5. 5.

Major mergers are not the dominant channel for triggering cluster formation in our model. The extremely high accretion rates onto dark matter halos at high-redshift are sufficient to trigger cluster formation without the aid of mergers (Fig. 9). At lower redshift, major mergers become increasingly important, especially for giant galaxies (Fig. 8). 6. 6.

The fraction of clusters formed during major mergers in our model is a factor of 4 lower than in recent simulations of galaxy formation (Li & Gnedin, 2019). This discrepancy may indicate a limitation of the analytical model in capturing detailed physics of GC formation. Nevertheless, both the analytical model and cosmological simulations predict similar distributions of GC formation epochs.

Acknowledgements

We thank Hui Li, Kareem El-Badry, Dan Weisz, Tom Zick, and Wren Suess for useful conversations and encouragement. We also thank the participants of the Lorentz center workshop “Formation of Stars and Massive clusters in Dwarf Galaxies over Cosmic Time” for useful discussions. This work made use of the matplotlib Python package (Hunter, 2007). This work was supported in part by the National Science Foundation through grant 1412144.