MODELING THE METALLICITY DISTRIBUTION OF GLOBULAR CLUSTERS

We follow the merging process of protogalactic dark matter halos using cosmological N-body simulations of three Milky Way-sized systems described in Kravtsov et al. (2004). The simulations were run with the adaptive refinement tree code (Kravtsov et al. 1997) in a 25 h⁻¹ Mpc box. Specifically, we use merger trees for three large host halos and their corresponding subhalo populations. The three host halos contain ∼10⁶ dark matter particles and have virial masses ∼10¹² M_☉ at z = 0. Two halos are neighbors, located 600 kpc from each other. The configuration of this pair resembles that of the Local Group. The third halo is isolated and is located 2 Mpc away from the pair. All three systems experience no major mergers at z < 1 and thus could host a disk galaxy like the Milky Way.

In addition to the host halos, the simulation volume contains a large number of dwarf halos that begin as isolated systems and then at some point accrete onto the host halo. Some of these satellites survive as self-gravitating systems until the present, while the rest are completely disrupted by tidal forces. We allow both the host and the satellite systems to form globular clusters in our model.

We adopt a simple hypothesis, motivated by the hydrodynamic simulation of Kravtsov & Gnedin (2005), that the mass in globular clusters, M_GC, that forms in a given protogalactic system is directly proportional to the mass of cold gas in the system, M_g. We define the corresponding mass fraction, f_g, of cold gas that will be available for star cluster formation in a halo of mass M_h as

$$\begin{equation} f_g \equiv {M_g \over f_b \, M_h}, \end{equation} \tag{ 1 }$$

where f_b ≈ 0.17 is the universal baryon fraction (Komatsu et al. 2010).

The gas fraction cannot exceed the total fraction of baryons accreted onto the halo, which is limited by external photoheating and depends on the cutoff mass M_c(z):

$$\begin{equation} f_{\rm {in}} = {1 \over (1 + M_c(z)/M_h)^{3}}. \end{equation} \tag{ 2 }$$

We use an updated version of the cutoff mass as a function of redshift (originally defined by Gnedin 2000) based on our fitting of the results of recent simulations by Hoeft et al. (2006), Crain et al. (2007), Tassis et al. (2008), and Okamoto et al. (2008):

$$\begin{equation} M_c(z) \approx 3.6\times 10^9 \, e^{-0.6\, (1+z)} \ h^{-1}\,{M}_{\odot }. \end{equation} \tag{ 3 }$$

Given the scatter in the simulation results and the numerical limitations of the modeling of gas physics, a reasonable uncertainty in this mass estimate is on the order of 50%. However, the resulting cluster mass and metallicity distributions are not very sensitive to the exact form of this equation. Note that Orban et al. (2008) provided an earlier revision of the equation for M_c(z); our current form is more accurate. If M_c(z) falls below the mass of a halo with the virial temperature of 10⁴ K, we set M_c(z) equal to that mass:

$$\begin{equation} M_{c, \rm min}(z) = M_4 \equiv 1.5 \times 10^{10} \Delta _{\rm vir}^{-1/2}\, {H_0 \over H(z)}\ h^{-1}\,{M}_{\odot }, \end{equation} \tag{ 4 }$$

where Δ_vir = 180 is the virial overdensity and H(z) is the Hubble parameter at redshift z. This criterion ensures that we only select halos that are able to cool efficiently via atomic hydrogen recombination lines.

Some of the baryons accreted onto a halo may be in a warm or hot phase (at T > 10⁴ K) and unavailable for star formation; thus, f_g < f_in < 1. We assume that only the gas in the cold phase (T ≪ 10⁴ K) is likely to be responsible for star cluster formation. The cold gas fraction f_g is calculated by combining several observed scaling relations. From the results of McGaugh (2005), the average gas-to-stellar mass ratio in nearby spiral and dwarf galaxies can be fitted as

$$\begin{equation} \frac{M_g}{M_*} \approx \left({\frac{M_*}{M_s(z)}}\right)^{-0.7}, \end{equation} \tag{ 5 }$$

where M_s is a characteristic scale mass, which we found to be M_s(z = 0) ≈ 4 × 10⁹ M_☉. This relation saturates at low stellar masses, where f_g cannot exceed f_in. At higher redshift the only information on the gas content of galaxies comes from the study by Erb et al. (2006) of Lyman break galaxies at z = 2. These authors estimated the cold gas mass by inverting the Kennicutt–Schmidt law and using the observed star formation rates. These estimates are fairly uncertain and model dependent. Within the uncertainties, their results can be fitted by the same formula but with a different scale mass: M_s(z = 2) ≈ 2 × 10¹⁰ M_☉. To extend this relation to all epochs, we employ a relation that interpolates the two values:

$$\begin{equation} M_s(z) \approx 10^{9.6 + 0.35 \, z}\,{M}_{\odot }. \end{equation} \tag{ 6 }$$

An additional scaling relation is needed to complement Equation (5) with a prescription for stellar mass as a function of halo mass. We compile it by combining the observed stellar mass–circular velocity correlation with the theoretical circular velocity–halo mass correlation. Woo et al. (2008) found that the stellar mass of the dwarf galaxies in the Local Group correlates with their circular velocities, which are taken as the rotation velocity for the irregular galaxies or the appropriately scaled velocity dispersion for the spheroidal galaxies. In the range 10⁷ M_☉ < M_* < 10¹⁰ M_☉, appropriate for the systems that may harbor globular clusters, the correlation is V_c ∝ M^{0.27 ± 0.01}_*. This can be inverted as M_* ≈ 1.6 × 10⁹ M_☉(V_c/100 km s⁻¹)^3.7.

Cosmological N-body simulations show that dark matter halos, both isolated halos and satellites of larger halos, exhibit a robust correlation between their mass and maximum circular velocity (e.g., Figure 6 of Kravtsov et al. 2004): V_max ≈ 100(M_h/1.2 × 10¹¹ M_☉)^0.3 km s⁻¹. This maximum circular velocity of dark matter is typically lower than the rotation velocity of galaxies because of the contribution of stars and gas. To connect the two velocities, we apply the correction $V_c = \sqrt{2} \, V_{\rm max}$, which reflects the observation that the mass in dark matter is approximately equal to the mass in stars over the portions of galaxies that contain the majority of stars. Then the equations in the last two paragraphs lead to M_* ≈ 5.5 × 10¹⁰(M_h/10¹² M_☉)^1.1 M_☉.

We also need to extend this local relation to other redshifts. Conroy & Wechsler (2009) matched the observed number densities of galaxies of given stellar mass with the predicted number densities of halos of given mass from z = 0 to z ∼ 2, averaged over the whole observable universe. They find that the stellar fraction f_* peaks at masses M_h ∼ 10¹² M_☉ and declines both at lower and higher halo masses. The range of masses of interest to us is below the peak, where we can approximate the f_* dependence on halo mass as a power law. The results from Figure 2 of Conroy & Wechsler (2009) are best fit by a steeper relation than we derived for the Local Group and also show significant variation with redshift at lower halo masses ∼10¹¹ M_☉: M_* ∝ M^1.5_h(1 + z)⁻² (there is much less variation with time around the peak at 10¹² M_☉, implying only a weak evolution in the total stellar density at z < 1). We adopt the same redshift dependence for our local relation, while using the shallower slope derived from Woo et al. (2008) because it deals with a population of halos in the mass range corresponding to the Milky Way progenitors. The stellar mass fraction of isolated halos is thus

$$\begin{equation} f_* = {M_* \over f_b \, M_h} \approx 0.32 \left({ \frac{M_h}{10^{12}\,{M}_{\odot }} }\right)^{0.1} (1+z)^{-2}. \end{equation} \tag{ 7 }$$

This relation steepens at low masses because of two additional limits on the gas and stellar fractions, which we impose to constrain the range of Equations (5)–(7) to the physical.

First, the sum of the gas and stars ("cold baryons") cannot exceed the total amount of accreted baryons in a halo:

$$\begin{equation} f_* + f_g \le f_{\rm {in}}. \end{equation} \tag{ 8 }$$

At each redshift, there is a transition mass M_h,cold, below which f_* + f_g = f_in and above which f_* + f_g < f_in. For masses M < M_h,cold (but not too low; see the next paragraph), we set f_g,revised = f_in − f_*, with f_* still given by Equation (7). We consider the baryons that are not included in f_g or f_* to be in the warm–hot diffuse phase of the interstellar medium.

Second, the ratio of stars to cold baryons, μ_* ≡ f_*/(f_* + f_g), is not allowed to increase with decreasing halo mass. For massive halos (M_h > M_h,cold), μ_* monotonically decreases with decreasing mass because of condition (5). At some intermediate masses M_h,μ < M_h < M_h,cold, μ_* continues to decrease but the gas fraction is reduced by condition (8). At M_h < M_h,μ, μ_* would reverse this trend and increase with decreasing halo mass because the cold gas is almost completely depleted. Such a reversal is unlikely to happen in real galaxies, which would not be able to convert most of their cold gas into stars. Therefore, for all masses M_h < M_h,μ we fix μ_* to be equal to the minimum value reached at M_h,μ. This affects both f_* and f_g.

We expect our stellar mass prescription to apply in the range of halo masses from 10⁹ to 10¹² M_☉, at least for the Local Group. However, this relation breaks when a halo becomes a satellite of a larger system. Satellite halos often have dark matter in the outer parts stripped by tidal forces of the host, while the stars remain intact in the inner parts. Unless the satellite is completely disrupted, we keep its stellar mass fixed at the value it had at the time of accretion, even though the halo mass may subsequently decrease.

The simultaneous effects of the above scaling relations are difficult to understand as equations. Figure 1 illustrates graphically the values of the gas and stellar fractions used in our model at various cosmic times. At z = 0 the gas fraction peaks for halos with M_h ∼ 3 × 10¹⁰ M_☉. At lower masses it is reduced by the amount of accreted baryons (Equation (8)), while at higher masses it is reduced by the gas-to-stars ratio (Equation (5)). The stellar fraction follows Equation (7) at high masses but drops faster at low masses because of constraint (8). For a Galaxy-mass halo, M_h ≈ 10¹² M_☉, our model gives M_* ≈ 5.5 × 10¹⁰ M_☉ and M_g ≈ 9 × 10⁹ M_☉. These numbers are consistent with the observed amount of the disk and bulge stars and the atomic and molecular gas in the Galaxy (from Binney & Tremaine 2008).

Zoom In Zoom Out Reset image size

At earlier epochs at all masses of interest, the gas fraction is higher and the stellar fraction is lower. There is a range of halos with M_h ≳ 10¹⁰ M_☉, which have an almost 100% gas fraction at redshifts z > 3. Such halos should be most efficient at forming massive star clusters.

We realize that our adopted relations for the evolution of the stellar and gas mass are not unique, as we are basing each fit on two data points. In order to test the sensitivity of our results to these assumptions, we consider alternative functional forms for these fits in Section 4.5. In particular, we give the stellar fraction a steeper dependence on halo mass and weaker dependence on cosmic time:

$$\begin{equation} f_{*,\rm alt} = 0.32 \left({\frac{M_h}{10^{12}\,{M}_{\odot }}}\right)^{0.5} (1+z)^{-1}. \end{equation} \tag{ 9 }$$

Such a slower evolution of the stellar mass is consistent with the observational studies of Borch et al. (2006), Bell et al. (2007), and Dahlen et al. (2007). The corresponding gas and stellar fractions are shown in Figure 2. Note that the amount of cold gas available for cluster formation is not strongly affected by this change (compare Figures 1 and 2).

Zoom In Zoom Out Reset image size

Having fixed the parameterization of the available cold gas, we then relate the gas mass of a protogalaxy to the combined mass of all globular clusters it can form within ∼10⁸ yr (the timescale of the simulation output). We based it on the rate derived in Kravtsov & Gnedin (2005):

$$\begin{equation} M_{\rm GC} = 3 \times 10^6 \,{M}_{\odot }\left({1 + p_2}\right) {M_g/f_b \over 10^{11} \,{M}_{\odot }}. \end{equation} \tag{ 10 }$$

An additional factor, 1 + p₂, allows us to boost the rate of cluster formation. Such a boost may be needed because we form new clusters only at arbitrarily chosen epochs corresponding to the simulation outputs. Unresolved mergers between the outputs may require p₂ > 0. In our model, we find the best fit to the Galactic metallicity distribution for p₂ ∼ 3 (see Table 1).

Note also that Equation (10) imposes the minimum mass of a halo capable of forming a globular cluster. Based on dynamical disruption arguments (Section 3) we track only clusters more massive than M_min = 10⁵ M_☉. Since we always have M_g < f_b M_h, in order to form even a single cluster with the minimum mass, the halo needs to be more massive than 10⁹ M_☉. For gas-rich systems at high redshift, M_GC ∼ 10⁻⁴ M_g/f_b ∼ 10⁻⁴ M_h.

Given the combined mass of all clusters to be formed in an event, M_GC, our procedure for assigning masses to individual clusters is as follows. We first draw the most massive cluster, which we call the nuclear star cluster, even though we do not have or use the information about its actual location within the host galaxy and it is not important for our current study. The mass assigned to the nuclear cluster, M_max, is derived from the assumed initial cluster mass function, dN/dM = M₀M⁻²:

$$\begin{equation} 1 = \int _{M_{\rm max}}^{\infty } \frac{{\it dN}}{{\it dM}} {\it dM}, \end{equation} \tag{ 11 }$$

which gives M_max = M₀. This normalization is constrained by the integral cluster mass:

$$\begin{equation} M_{\rm GC} = \int _{M_{\rm min}}^{M_{ \rm max}} M \frac{{\it dN}}{{\it dM}} {\it dM} = M_{\rm max} \ln {M_{\rm max} \over M_{\rm min}}. \end{equation} \tag{ 12 }$$

The power-law initial mass function (IMF) agrees both with the observations of young star clusters and the hydrodynamic simulations. After the nuclear cluster is drawn, the masses of smaller clusters are selected by drawing a random number 0 < r < 1 and inverting the cumulative distribution: r = N(<M)/N(<M_max):

$$\begin{equation} M = \frac{M_{\rm min}}{1 - r \left({1 - M_{\rm min}/M_{\rm max}}\right)}. \end{equation} \tag{ 13 }$$

We continue generating clusters until the sum of their masses reaches M_GC.

The formation of clusters is triggered by a gas-rich major merger of galaxies, which includes mergers of satellite halos onto the main halo as well as satellite–satellite mergers. New clusters form when the halo mass at the ith simulation output exceeds the mass at the previous output by a certain factor, and the cold gas fraction exceeds a threshold value at the same time:

$$\begin{equation} {\tt case\hbox{-}1:} \quad M_{h,i} > (1+p_3)\ M_{h,i-1} \quad {\rm and} \quad f_g > p_4. \end{equation} \tag{ 14 }$$

Also, we require that the maximum circular velocity not decrease in this time step to ensure that the mass increase was real rather than a problem with halo identification. We have experimented with a more relaxed criterion for the main halo than for satellite halos, with p_3,main < p_3,sat, but did not find a significantly better fit to the mass or metallicity distributions. We therefore keep a single value of p₃ for all halos.

For some model realizations, we consider an optional alternative channel for cluster formation without a detected merger, if the cold gas fraction is very high:

$$\begin{equation} {\tt case\hbox{-}2:} \quad f_g > p_5, \end{equation} \tag{ 15 }$$

where the threshold p₅ is expected to be close to 100%. This channel allows continuous cluster formation at high redshift when the galaxies are extremely gas-rich. High-redshift galaxies are probably in a continuous state of major and minor merging, but because of their lower masses it is more difficult to detect such mergers in the simulation. Additional motivation for this channel follows from some nearby starburst galaxies that are forming young massive clusters despite appearing isolated. Case-2 formation is allowed only for isolated halos before they are accreted into larger systems and become satellites. The epoch of accretion is defined by the last timestep before the orbit of the subhalo falls permanently within the virial radius of its host.

Our model sample combines clusters formed in the main halo and in its satellites, either surviving or disrupted. We exclude clusters from the satellites that have a galactocentric distance at z = 0 greater than 150 kpc, which is the largest distance of a Galactic globular cluster. We apply the criteria for cluster formation at every timestep of the simulation (every ∼10⁸ yr) for each of the three main halos and their satellite populations. The rate of cluster formation per every merger event is therefore approximately M_GC/10⁸ yr.

In order to compare the distribution of clusters obtained from our analysis to the distribution of Galactic globular clusters, we normalize the total number of model clusters by the ratio of the Galaxy mass to the simulated halo masses at z = 0:

$$\begin{equation} N_{\rm normalized} = N_{\rm model}\ \frac{M_{\rm {MW}}}{M_{h1} + M_{h2} + M_{h3}}. \end{equation} \tag{ 16 }$$

We take M_MW = 10¹² M_☉ and use M_h1 = 2.37 × 10¹² M_☉, M_h2 = 1.77 × 10¹² M_☉, and M_h3 = 1.70 × 10¹² M_☉ from Kravtsov et al. (2004).

The iron abundance is assigned to each model cluster according to the estimated average metallicity of its host galaxy. We obtained the latter from the mass–metallicity relation for dwarf galaxies of the Local Group at z = 0 as formulated by Woo et al. (2008):

$$\begin{equation} \mathrm{[Fe/H]}_0 = -1.8 + 0.4\log {\left({\frac{M_*}{10^6 \,{M}_{\odot }}}\right)}. \end{equation} \tag{ 17 }$$

In fact, the same fit is valid for the smallest, ultrafaint dwarfs studied by Kirby et al. (2008). Thus, we apply this relation to all protogalactic systems in our simulation volume.

We also include the evolution of this relation with cosmic time, based on the available observations of Lyman-break galaxies at z ≈ 2 (Erb et al. 2006), Gemini Deep Survey galaxies at z ≈ 1 (Savaglio et al. 2005), and cosmological hydrodynamic simulations that provide the average metallicity of galaxies (Brooks et al. 2007; Davé et al. 2007):

$$\begin{equation} \mathrm{[Fe/H]}(t) \approx \mathrm{[Fe/H]}_0 - 0.03 \left({t_0 - t \over 10^9 \ {\rm yr}}\right). \end{equation} \tag{ 18 }$$

While this temporal evolution is probably real, it can change the metallicity for the same stellar mass by at most 0.36 dex in 12 Gyr. This amount is smaller than the 0.4 dex change of [Fe/H] due to each factor of 10 variation of stellar mass. In our model, globular clusters form in protogalaxies with a range of stellar masses of several orders of magnitude (see Figure 9 below).

The observed mass–metallicity relation for a large sample of galaxies observed by the Sloan Digital Sky Survey has an intrinsic scatter of at least 0.1 dex (e.g., Tremonti et al. 2004). We account for it, as well as for possible observational errors, by adding a Gaussian scatter to our calculated [Fe/H] abundances with a standard deviation of σ_met = 0.1 dex. The exact value of this dispersion is not important and can go up to 0.2 dex without affecting the results significantly.

Using Equations (17) and (18) along with the procedures of Section 2.2, we can generate a population of star clusters with the corresponding masses and metallicities. The model contains two random factors: the scatter of metallicity and the individual cluster masses assigned via Equation (13). We sample these random factors by creating 11 realizations of the model with different random seeds. Each realization combines clusters in all three main halos. Taking into account that the halos are about twice as massive as the Milky Way, the expected number of clusters in each model realization is ∼150(the observed number) × (2.37 + 1.77 + 1.7) ≈ 870. The total set of all 11 realizations includes ∼9500 clusters. For the purpose of conducting statistical tests on the distributions of cluster mass and metallicity, we consider each realization separately and then take the median value of the calculated statistic.

For convenience, we provide a list of the most important equations we used in the model in Table 2.

Overall, our model has five adjustable parameters (Table 1). To explore possible degeneracies among these parameters, and to find the parameter set that produces the best-fitting metallicity distribution, we set up a grid of models in which each of the parameters was varied within a finite range of values. The range was taken to be large enough to explore all physically relevant values of each parameter.

The boost for cluster formation, p₂, varied from 0 to 5. For consistency with the rate derived in the hydrodynamic simulation of Kravtsov & Gnedin (2005), we aimed to keep this parameter at low values.

The minimum mass ratio for mergers, p₃, varied between 0.15 and 0.5. It is consistent with typical major merger criteria used in the literature (e.g., Beasley et al. 2002 use p₃ = 0.3).

The cold gas fraction required for cluster formation during a merger, p₄, could be relatively low but non-zero, so that we considered 0 < p₄ < 0.2. This threshold parameter accounts for why disk galaxies like the Milky Way are still forming stars despite a low gas fraction, while ellipticals are not.

The gas fraction for case-2, p₅, has to be very high—above 90%, as our prescription predicts that many halos have a very high gas fraction at high redshift and could overproduce blue globular clusters (as was the case in the Beasley et al. 2002 model).

We considered several values for σ_met but found that a value of 0.2 or higher smeared out the peaks in the metallicity distribution, while a value of 0 failed to fill the extreme ends of the distribution. We therefore include only three values in our search, σ_met = 0, 0.1, 0.2.

We find the best-fit model by searching through the multi-parameter space and maximizing the K-S probabilities of the metallicity distribution, P_KS,Z, and the mass function, P_KS,M, being consistent with observations. The likelihood function also contains additional factors that force the parameters toward the values that we consider ideal. We require the model to produce the observed number of clusters, N ≈ 150, scaled by the host galaxy mass as in Equation (16). We wish to maximize the fraction of clusters formed in the main disk, f_disk, to be consistent with the observed spatial distribution (Section 5). We penalize the likelihood function for large values of p₂ and for any young clusters formed after t = 10 Gyr, N_after10. We also wish to minimize the fraction of clusters formed through the case-2 channel, f_case2, for simplicity of the model. Finally, we want to increase the likelihood of the metallicity distribution being bimodal, as characterized by the dip test, P_dip, which we discuss later in Section 4.4. The actual likelihood function that we maximize is given by

$$\begin{eqnarray} \log {\cal L} &=& \log {P_{\rm {KS},Z}} + 0.3\log {P_{\rm {KS},M}} -[(N-150)/30]^2 \nonumber \\ && +\log {f_{\rm disk}} - 0.15 p_2 - 20 N_{\rm after10}/N - 0.4 f_{\rm case2}\nonumber \\ && + 3 \log {P_{\rm dip}}. \end{eqnarray} \tag{ 26 }$$

The coefficients for each term were adjusted heuristically until we found that their relative weights matched our expectation for selecting acceptable distributions. The "best-fit" distribution that maximizes ${\cal L}$ is therefore a subjective fiducial model that we use to illustrate how the bimodality may arise. We then look at how many model realizations are similar to the "best-fit" for other possible values of the parameters.

Figure 6 shows the predicted best-fit metallicity distribution of model clusters and the observed distribution of Galactic globular clusters, both metal-poor and metal-rich. Note that we require our model to have the same formation criteria for both cluster populations; we do not explicitly differentiate between the two modes. The only variable is the gradually changing amount of cold gas available for star formation, yet the model predicts two peaks of the metallicity distribution, centered on [Fe/H] = −1.54 and [Fe/H] = −0.58, in remarkable agreement with the observations. The standard deviation of the red peak is 0.24 dex and that of the blue peak is 0.32 dex.

Zoom In Zoom Out Reset image size

The K-S probability test of the model and data samples being drawn from the same distribution is P_KS,Z = 80%, that is, they are fully consistent with each other. The number of surviving clusters is N = 147, also matching the observations. Even though our current model is extremely simple, this bimodality is reproduced naturally, without explicit assumptions about truncation of the production of metal-poor clusters at some early epoch or about the formation of metal-rich clusters in a merger of two spiral galaxies.

We find that the main halo contributes more significantly to the red peak than it does to the blue peak (Figure 7). In particular, clusters with the highest [Fe/H] appear to have been formed primarily by late merging into the main halo.

Zoom In Zoom Out Reset image size

The fraction of clusters formed via the case-2 channel is f_case2 = 22%. These clusters produce a single-peaked distribution of blue clusters. In contrast, clusters formed in major mergers contribute to both red and blue modes, in about equal proportions. We return to this point in the discussion of globular cluster systems of elliptical galaxies in Section 7.

Clusters that formed after z = 2 constitute the bulk of the red peak and contribute little to the blue peak in the metallicity distribution (Figure 8). The strength of this result implies that the gas reservoir and the rate of hierarchical merging at intermediate redshifts are conducive to the creation of red clusters. This result lends itself well to the idea that the simulation of Kravtsov & Gnedin (2005) was only able to reproduce the metal-poor population of globular clusters because the simulation was stopped at z ≈ 3.

Zoom In Zoom Out Reset image size

Our prescription links cluster metallicity to the average galaxy metallicity in a one-to-one relation, albeit with random scatter. Since the average galaxy metallicity grows monotonically with time, clusters forming later have on the average higher metallicity. The model thus encodes an age–metallicity relation, in the sense that metal-rich clusters are younger by several Gyr than their metal-poor counterparts. This relation is required in the model to reproduce the observed metallicity distribution because very old galaxies cannot produce high enough metallicities. However, Figure 8 shows that clusters of the same age may differ in metallicity by as much as a factor of 10, as they formed in the progenitors of different mass.

Available observations of the Galactic globular clusters do not show a clear age–metallicity relation, but instead indicate an age spread increasing with metallicity (De Angeli et al. 2005; Marín-Franch et al. 2009; Dotter et al. 2010; Forbes & Bridges 2010). Red clusters have a younger mean age overall and may be as young as τ ≈ 7 Gyr. Our model does not appear to be in an obvious conflict with this trend. We define cluster age as τ ≡ t₀ − t_f, where t_f is the time of formation. We find the mean age of 11.7 Gyr for the blue population and 6.4 Gyr for the red population, with a standard deviation of 1.3 Gyr and 2.7 Gyr, respectively. More accurate dating of the Galactic and extragalactic clusters is needed to falsify the predicted age–metallicity trend.

Distributions of the cluster formation time and environment in the fiducial model are shown in Figure 9. The age distribution, which peaks strongly between 11 and 13 Gyr, demonstrates that the majority of our clusters are still very old and falls in line with the observed perception of globular clusters. However, the distribution of formation redshift appears remarkably flat in the range 1 < z < 7, emphasizing that the clusters were not formed in a single event but rather through the continuous process of galaxy formation. Few clusters were formed prior to the era of reionization, as sufficiently large quantities of gas could not be condensed to meet the mass threshold for cluster formation at redshifts z > 9. The distributions of the total and stellar mass of the host galaxies extend over 3 orders of magnitude. Their extended high-mass tails contribute to the strength of the red peak, as the most massive halos would form the most metal-rich clusters.

Zoom In Zoom Out Reset image size

Globular clusters form much earlier than the majority of field stars. Figure 10 shows the fraction of galaxy stellar mass locked in massive star clusters, normalized for convenience as 10³ M_GC/M_*. To calculate this ratio, we summed over all protogalactic systems that would end up within 150 kpc of the galaxy center at z = 0, regardless of their location at earlier times. Thus, it represents a global cluster formation efficiency in a Milky Way-sized environment. Specific realizations of the model differ in detail in the three host halos by as much as a factor of 2. This scatter is shown by the shaded region on the plot. The globular cluster mass includes their continuous formation and the mass loss due to the dynamical evolution. A striking prediction of the model is a very high cluster fraction at early times, near t = 1 Gyr, of M_GC/M_* ≈ 10%–20%. Star cluster production may have been a dominant component of galactic star formation at z > 3. By t = 3 Gyr (z ≈ 2), the cluster fraction drops to only a few percent, as expected for a galaxy undergoing active star formation. At the current epoch, massive star clusters make up less than 0.1% of the stellar mass. The predicted ratio is progressively more uncertain at higher redshifts because it relies on our extrapolated prescription for the galactic stellar mass. The low-redshift prediction should be robust. We also show a variant of the specific frequency parameter related to the number of clusters, T ≡ N/(M_*/10⁹ M_☉), introduced by Zepf & Ashman (1993). It shows a similar decline with time, reaching T ≈ 2 at the present.

Zoom In Zoom Out Reset image size

These global cluster formation efficiencies agree with many observations across galaxy types. Rhode et al. (2005) find T ∼ 1 for both red and blue clusters in the field and group spiral galaxies. This parameter increases with the galaxy mass. In the Virgo cluster, Peng et al. (2008) find T ∼ 5 for galaxies in the mass range appropriate for the Milky Way. McLaughlin (1999) estimated the cluster mass fraction in both spiral and elliptical galaxies to be M_GC/(M_* + M_g) ≈ 0.0026 ± 0.0005. This is larger than what we find by a factor of several, but we count in M_* all stars out to 150 kpc, which includes some satellite galaxies as well as the host. Therefore, both predicted cluster efficiencies at z = 0 are reasonable. Their rise at high redshift is an interesting prediction of the model.

The model also shows that the globular cluster system overall is more metal-poor than the stars in disrupted satellites, which are expected to form a stellar spheroid of the Galaxy. We calculated the mass-weighted metallicity of stars formed in the disrupted satellites of all three main halos (using Equations (7), (8), (17), (18)). This calculation bears all the uncertainty of our extrapolated time evolution of the stellar fraction and mass–metallicity relation, but nevertheless provides a useful estimate. We find the tail of halo star metallicities as low as the most metal-poor globular clusters, but the overall stellar distribution peaks around [Fe/H] ≈ −0.3. A very similar situation is observed in NGC 5128 and discussed by Harris (2010). In our model, the majority of globular clusters form before the bulk of the field stars and therefore acquire lower metallicities. For comparison, the metallicity of stars in surviving satellite galaxies peaks around [Fe/H] ≈ −0.8 and forms an intermediate population between the clusters and the field.

Despite our attempts to incorporate it as a major penalty in the likelihood statistic, we were unable to completely eliminate the phenomenon of young massive star clusters. Interestingly, these clusters did not originate in the main galactic disks. All clusters younger than 5 Gyr formed in satellite halos in the mass range ∼10¹⁰–10¹¹ M_☉, at distances 40–100 kpc from the center. Although the proper sample of the Galactic globular clusters does not contain any young clusters, there are several young massive clusters in M31 whose ages were confirmed both from the visual and UV colors (Fusi Pecci et al. 2005; Rey et al. 2007) and from the integrated-light spectroscopy (Puzia et al. 2005). The actual analogs of young model clusters may be found in the Large Magellanic Cloud (LMC), which hosts globular clusters with a wide range of ages and continues to form clusters unti now. There may even exist young star clusters with masses ∼10⁵ M_☉ in the Galactic disk, hidden behind tens of visual magnitudes of extinction but revealing themselves through free–free emission of their ionization bubbles (Murray & Rahman 2010). Massive star cluster formation at late times thus paints a picture consistent with the idea that today's super star clusters are destined to become observationally equivalent to globular clusters, as envisioned by Ashman & Zepf (1992) and Harris & Pudritz (1994).

A separate criterion for the formation of clusters in extremely gas-rich systems (case-2) is not necessary for achieving a good fit to the observed metallicity distribution. Though we feel that the inclusion of the case-2 formation channel in the model is both useful and physically motivated, it takes away from the elegance of using only resolved mergers as the lone formation mechanism. It turns out that the main benefit of allowing clusters to form via case-2 is seen in the mass function of surviving clusters. The high-mass end of the mass function matches the observations better if massive halos (primarily the main halo) are allowed to form as many clusters as possible at early times.

We searched the model grid without the case-2 channel formation, by setting p₅ = 1, and found an almost equally good metallicity distribution as in the fiducial model. Figure 11 shows that this distribution also appears bimodal and completely consistent with the data. The K-S probability is P_KS,Z = 92%. In fact, even the mass function is only marginally less consistent, P_KS,M = 2.0% versus 7.4% in the fiducial model. The parameters used to obtain this distribution were p₂ = 2.85, p₃ = 0.16, p₄ = 0.04, p₅ = 1, σ_met = 0.1.

Zoom In Zoom Out Reset image size

The fiducial distribution discussed above is not unique among our results in its ability to match the observations. Significant degeneracy exists among combinations of the model parameters that produce metallicity distributions consistent with the Galactic sample. Many models within the grid have sufficiently high K-S probabilities. In this section, we explore which regions of the parameter space produce models similar to our best fit.

First, let us motivate the use of the likelihood function given by Equation (26) as opposed to using a standard statistical test to select the best fit. In the early stages of development of our model, we relied on the K-S test alone to help us understand the range of parameters that produce metallicity distributions that match the observations. However, once the model was completed, it became apparent that the K-S test alone is not powerful enough for the analysis of the results. This is clearly demonstrated in Figure 12, which shows the value of the K-S probability P_KS,Z as a function of p₂ and p₃ across their respective ranges in the grid. Each point represents the maximum possible P_KS,Z for the given values of p₂ and p₃ with the other parameters free to vary within the grid. This is done to best represent the full extent of the five-dimensional parameter space within a two-dimensional slice. Statistically, any distribution with P_KS,Z > 10% cannot be ruled out with confidence, implying that almost the entire range of our parameters can produce statistically consistent distributions. In addition, although some regions of the parameter space have higher values of P_KS,Z than others, there is no clear pattern in the contours to help us understand the required physics of star cluster formation within our semi-analytical recipe.

Zoom In Zoom Out Reset image size

In comparison, Figure 13 shows contours of the value of the likelihood function from Equation (26), using the same scheme described above to maximize the value at each point. The shape of these contours demonstrates that p₂ and p₃ are degenerate in their ability to produce good distributions. The degeneracy can most easily be understood by noting that these parameters directly affect the total number of clusters: p₂ controls the cluster formation rate per merger, while p₃ selects eligible mergers. It is therefore expected that the contours show a correlation at high levels of the likelihood function, as the statistic depends sensitively on the total number of clusters.

Zoom In Zoom Out Reset image size

Figure 14 shows the same type of contours as Figure 13, but only for distributions with p₅ = 1. Not using the case-2 channel reduces the range of the parameter space where good distributions are found. In particular, compared to the previous plot, Figure 14 lacks any viable models with p₃ > 0.3. Given the tight and steep correlation in this plot, it is likely that larger values of p₃ would require very high p₂ > 5, which may violate current observational constraints on the cluster formation efficiency. However, Figure 11 demonstrates that a good model without the case-2 channel can still be found with reasonably small values of p₂ and p₃.

Zoom In Zoom Out Reset image size

To understand the sensitivity of the likelihood function to individual parameters, we also considered one-dimensional slices of the parameter space around the fiducial model, this time allowing only one parameter to vary at a time. Figure 15 illustrates how the sharp peaks of ${\cal L}$ allow us to select the best model more accurately than on the basis of P_KS,Z alone. Particularly, as a function of p₂ and p₃, P_KS,Z varies slowly over the entire range of the grid. On the other hand, p₄ and p₅ must stay within a small range of their fiducial values in order to achieve acceptable values of either P_KS,Z or ${\cal L}$.

Zoom In Zoom Out Reset image size

Figure 16 shows variation of the metallicity distribution when individual parameters deviate from their fiducial values. Each parameter can change the shape of the metallicity function and the number of clusters. The effects of varying p₂ and p₃ are almost opposite, reflecting the degeneracy in the likelihood contours. In particular, a smaller p₃ accommodates more minor mergers, which allow massive hosts to form more metal-rich clusters as well as some metal-poor clusters. Decreasing p₅ allows more clusters to form through the case-2 channel; most of such clusters are metal-poor. The major role of p₄ appears to govern the extent of the most metal-rich clusters—lower threshold gas fraction allows clusters to form in the later, more enriched environments of massive hosts.

Zoom In Zoom Out Reset image size

Figure 17 illustrates the response of the metallicity distribution to simultaneous variations of model parameters. First, we plot two distributions where we changed p₃ to 0.15 and 0.3 while keeping the other parameters fixed. The width of the metal-poor peak broadens as p₃ is lowered, indicating that a wider range of halos in the early universe were able to produce clusters. Raising p₃ has the opposite effect. Note that the locations of the two peaks are remarkably robust to these changes. Staying at p₃ = 0.15, we set p₅ = 1 to eliminate the case-2 channel and set p₄ = 0 to allow even gas-poor massive halos at low redshift to form clusters. The result (long-dashed line) is a distribution with a much broader metal-rich peak, which extends well past the maximum metallicity of the fiducial model. The dot-dashed line represents a corresponding change to p₄ = 0.08 and p₅ = 0.96 for the p₃ = 0.3 model. In this case, the metal-rich peak is severely depleted and remains only as an extended tail of a single-peaked, metal-poor distribution dominated by case-2 clusters. These distributions are just some of the realizations of our model that were rejected due to their low values of ${\cal L}$. All of them have features that conflict with the observed metallicity distribution in the Galaxy.

Zoom In Zoom Out Reset image size

The K-S statistic measures the overall consistency of the model and observed metallicity distributions, but not specifically bimodality or multimodality within the distributions. In order to address the particular issue of modality, we employ two additional statistical tests, described in the Appendix.

The Gaussian mixture modeling (GMM) test indicates that the fiducial distribution is bimodal at a high level of significance (better than 0.1%). The peak metallicities of both modes and their widths are close to the observed values and agree with them within the errors. Both samples easily appear bimodal to the eye because the modes are well separated, with the dimensionless peak separation ratio D > 3 (see Equation (A3)). However, as we discuss in the Appendix, the GMM test is sensitive to the assumption of Gaussian modes. It may indicate a highly statistically significant split into two modes when the distribution is truly unimodal but skewed. For faster and more robust model selection we consider another test of multimodality.

The dip test compares the cumulative input distribution with the best-fitting unimodal distribution. The maximum distance between the two corresponds to a dip in the differential distribution. The dip test for the observed Galactic clusters indicates that the distribution is 90% likely to be not unimodal. When applied to our fiducial model, the dip test implies it is 99% likely to be not unimodal. However, there is a caveat that the probability of the dip test depends on the number of objects in the sample, similar to the K-S test. The higher significance of the model result does not mean that the model is actually more bimodal than the data because we used all 11 random realizations of the model as a combined sample to evaluate the dip test. While this is not a fair comparison to the data, it allows us to differentiate efficiently among alternative models.

We ran the dip test for all models on the grid in a manner similar to the likelihood statistic. The most interesting result of the dip statistic comes from one-dimensional slices of the parameter space. Considering only models with the normalized number of clusters in the range 140 < N < 160, we binned the distributions according to the values of the four parameters and found the median and quartiles of P_dip in each bin. Figure 18 shows several trends. (1) Distributions with low formation rate p₂ are unlikely to be bimodal. The 75th percentile of P_dip increases systematically with p₂ in the range 2 < p₂ < 3, but plateaus for p₂ > 3. (2) The most-bimodal distributions require p₃ to be small enough to allow for merger ratios 1:5 to trigger cluster formation. Between p₃ = 0.2 and p₃ = 0.5, the lower p₃ is, the better. However, mass ratios lower than 1:6 may dilute bimodality. (3) The gas fraction threshold p₄ should be under 10% for ideal bimodality, to include mergers of massive galaxies. (4) The fraction p₅ has to be close to 1, implying that case-2 negatively affects bimodality. A conclusion from this plot is that bimodality appears in a significant number of model realizations, for a wide range of parameters. At the same time, a similarly large number of realizations are unimodal.

Zoom In Zoom Out Reset image size

The metallicity distribution is bimodal if metal-rich clusters constitute a significant subset of all clusters. Thus, the fraction of red clusters, f_red, is a simple proxy for bimodality. Indeed, we find a strong correlation between f_red and P_dip. The red fraction follows similar, but weaker, trends with the model parameters to those shown in Figure 18. The median red fraction correlates most strongly with p₅, increasing from f_red = 16% for p₅ = 0.96 to f_red = 32% for p₅ = 1. The red peak is significantly stronger without case-2 clusters.

We note that the dip test, unlike the K-S test, does not depend on comparing the model distribution to the Galactic sample. Therefore, the trends for bimodality derived from the dip test should apply to other globular cluster systems. We anticipate that bimodality would likely arise if we applied our model with the parametric constraints stated above to any cosmological N-body simulation that follows the mass assembly history of a large galaxy. Further discussion of applying our model in different galactic environments follows in Section 7.

In order to investigate the underlying cause of bimodality, we examined various properties of merger events. A merger event is defined as any time in a halo's track when it meets the criteria for case-1 formation. An important requirement here is the minimum mass of cold gas needed to produce a cluster that would survive dynamical disruption. Through Equation (10), we see that a cluster mass M > 2 × 10⁵ M_☉ requires M_g > 3 × 10⁸ M_☉. This constraint significantly reduces the number of eligible mergers. We considered the distributions of halo mass, lookback time, and metallicity (without additional dispersion) for all relevant merger events. We find that relatively few mergers happen in the space of high metallicity, high mass, and late time. Almost half of the mergers (44%) take place before τ = 12 Gyr, and only 24% of the mergers happen in the last 10 Gyr. If we also counted the events that led to now-disrupted clusters, these numbers would spread even further to 53% and 17%, respectively. Nevertheless, the recent mergers stand out for two reasons: each such event creates more clusters, and these younger clusters have a better chance of surviving the dynamical disruption than the older clusters. Since the number of clusters formed in each merger is positively correlated with the galaxy mass, the few stochastic super-massive mergers with high metallicity are likely to produce a significant number of clusters, which would separate the red peak from the blue peak.

We also considered that cluster bimodality may be linked to the mass ratios in the merger events. "Major" and "minor" mergers have been proposed to play different roles in galaxy formation, so it is conceivable that different types of star clusters may be formed depending on the merger ratio. In Figure 19, we plot cumulative metallicity distributions for clusters grouped according to the mass ratios in the cluster-forming merger event. Mergers with the closest masses, ΔM_h/M_h > 0.5, contribute 48% of case-1 clusters, while the two lower mass ranges each contribute equal portions of the rest. Running a K-S test on these distributions revealed that they formally represent statistically different populations. However, there is no clear-cut range of metallicities where one type of merging is exclusively producing all of the clusters, and the overall shapes of the distributions are similar. This uniformity suggests that bimodality is a natural consequence of hierarchical cluster formation regardless of the exact definition of a "major" merger.

Zoom In Zoom Out Reset image size

Figure 20 shows the number of models from the grid that fall into particular ranges of the dip probability and the ratio of case-2 clusters to case-1 clusters, N₂/N₁. The models are restricted to have the normalized number of clusters 140 < N₁ + N₂ < 160. The region with the highest density of models is in the lower-right corner of the plot, corresponding to high P_dip and low N₂/N₁. Low values of P_dip are not significant since they cannot reject a unimodal distribution. Effectively, bimodality requires N₂/N₁ ≲ 0.5. At the significance level of P_dip = 90%, corresponding to the observed distribution, 38% of the grid models are bimodal if N₂/N₁ < 0.3. This fraction drops to only 15% for 0.3 < N₂/N₁ < 1, and then further to 9% for N₂/N₁ > 1. These statistics confirm that bimodal populations appear only when the case-2 channel is a secondary formation mechanism.

Zoom In Zoom Out Reset image size

Another part of the explanation of bimodality of the surviving clusters is due to the dynamical evolution. Most of the disrupted clusters were old and blue. If we add these disrupted clusters to the metallicity histogram in the fiducial model, the blue peak rises by a factor of 2. The red peak remains virtually unaffected, since the more recently formed red clusters are less subjected to dynamical disruption. We ran the dip test on all grid distributions, including both surviving and disrupted clusters. Among the models with the number of surviving clusters in the range 100 < N < 200, few distributions have P_dip > 50% and none has P_dip > 80%. This means that there is virtually no bimodality. Indeed, the distributions appear almost entirely unimodal, with the peak in the blue end and nothing more than a tail in the red end. This leads to a prediction that late-type galaxies, which have more continuous cluster formation than early-type galaxies, may be less likely to exhibit bimodal cluster populations.

As alluded to in Section 2, some equations that we used in the prescriptions for the stellar mass and the cold gas fraction were based on only a few observed points. Currently there is limited observational or theoretical understanding of how these functions should behave at high redshift, which is the period of primary interest for our study. Below we consider some alternatives for these prescriptions.

The stellar fraction that we adopted from Woo et al. (2008) is well motivated by Milky Way dwarf galaxies at the present epoch, but the abundance-matching models such as the one by Conroy & Wechsler (2009) predict a steeper dependence on halo mass in the range 10⁸ M_☉ < M_* < 10¹⁰ M_☉. Additionally, the redshift dependence of this relation is uncertain, and recent observational surveys (Borch et al. 2006; Bell et al. 2007; Dahlen et al. 2007) have advocated a slower evolution than the (1 + z)⁻² adopted in our model. To accommodate this uncertainty, we re-ran the entire parameter grid with Equation (9) instead of Equation (7). The corresponding contour plot is shown in Figure 21 and the best-fit metallicity distribution is shown in Figure 22. This best fit is capable of reproducing the observed metallicity (P_KS,Z = 49%) and mass distributions (P_KS,M = 9.5%), similar to our fiducial model. The acceptable range of the parameter space is narrower and shifted toward higher values of p₂, as the steeper mass slope otherwise prevents low-mass halos from forming a sufficient number of clusters. Nevertheless, this alternative prescription still leads to a significant chance of a bimodal metallicity distribution.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Current observational constraints on the gas fraction at high redshift are even more uncertain. We adjusted the fit by altering the scale mass in Equation (6). As alternatives to the fiducial model, we considered a redshift power law, M_s = M_s,0(1 + z)², and an inverse time dependence, M_s = M_s,0(t/t₀)⁻¹. Both of these relations resulted in lower gas fractions during the high-redshift epoch when most globular clusters should form. The gas fractions were too low for case-2 formation for any reasonable choice of p₅. More importantly, the low gas masses did not allow even the most massive halos to form star clusters until intermediate redshifts. Therefore, none of these fits is a viable alternative to the fiducial model.

The same problems appeared if we took the simplistic approach of holding the gas fraction constant for all halo masses at all times. This idea was initially considered to see if we could generate simple results based only on halo merger histories without speculating on the baryonic physics. We quickly realized that this approach was not going to work. Setting f_g too low effectively prevents cluster formation at high redshift, when blue globular clusters are expected to form, as most halos cannot build up sufficient mass to overcome the minimum mass required to form a single massive star cluster (discussed in Section 2.2). A constant gas fraction that is too high, on the other hand, presents obvious unphysical predictions at low redshift, and in particular would drastically overpredict the number of young clusters, forcing us to arbitrarily truncate their formation.

We also considered an alternative parameterization of the gas fraction, suggested by Stewart et al. (2009). They took the same observational constraints as us, but fitted them as

$$\begin{equation} \frac{M_g}{M_*} = 0.04 \, \left({M_* \over 4.5 \times 10^{11}\,{M}_{\odot }}\right)^{-0.59 (1+z)^{0.45}}. \end{equation} \tag{ 27 }$$

This formula predicts so much cold gas at high redshift that many low-mass halos would be able to form clusters via the case-2 channel for any p₅ < 1. If we completely disable case-2 formation and use the above prescription for the gas mass, we find many model realizations consistent with the observed metallicity distribution. This prescription differs from our fiducial choice in that it produces considerably more young clusters and achieves less clear metallicity bimodality. The maximum value of the likelihood function attainable with this prescription is approximately half of the value for the fiducial prescription. Nevertheless, it could still produce acceptable globular cluster results.

In addition to changing the formulation of the fits, we investigated the effect of adding a random Gaussian dispersion with standard deviation σ_fits to the right-hand sides of Equations (3), (5), and (7) to reflect their intrinsic scatter as well as observational uncertainty. Different random values for the three scatters are generated for each halo at each timestep, but we always force the condition given by Equation (8) on the total baryon content. For simplicity, we used the same magnitude of σ_fits for the scatter added to all three equations simultaneously. We re-ran the parameter search grid using σ_fits = 0.1, 0.2, and 0.3 dex. For each value of σ_fits, we were still able to find models with high values of P_KS,Z and overall likelihood statistic, although these values decline with the increasing amount of scatter. The metallicity probability varies from 49% to 16% to 6%, for σ_fits = 0.1, 0.2, 0.3 dex, respectively.

As an alternative to scatter in the cutoff mass (Equation (3)) with a fixed functional form for the gas fraction (Equation (2)), we tried adding scatter to Equation (2) while keeping Equation (3) fixed. Adding scatter to f_in allows the gas fraction to exceed the threshold p₅ much more easily and to produce too many case-2 clusters. To avoid unphysical results, we analyze only the results for the case-1 formation channel. In this case, we find the best-fit models with P_KS,Z = 47%, 13%, and 5%, for σ_fits = 0.1, 0.2, 0.3 dex, respectively. These models are still consistent with the observed Galactic distribution.

The addition of scatter as described above has two systematic effects on any individual realization of the metallicity distribution: the high metallicity tail is extended even further and the height of the blue peak is damped relative to the case with no scattering. The former effect is due to the possibility of drawing higher values of M_* and hence higher [Fe/H]. The latter effect arises from the enforcement of Equation (8), which prevents gas-rich halos at high-redshift from gaining any extra gas from the positive scatter in Equation (5); on the other hand, negative scatter can prevent some of these halos from being eligible for case-2 formation. Accordingly, the best distributions with higher values of σ_fits were found for models with low values of p₅. We note that the dip probability in most realizations is not strongly affected by the new scatter, implying that the actual smearing of the peaks is not significant and bimodality is preserved.

In Section 3, we noted that the expression for the evaporation rate (Equation (23)) contains some inherent parameters. Here, we explore alternative disruption models with different values of ξ_e and δ within the fiducial formation prescription.

The effect of decreasing ξ_e is simply to reduce the number of clusters that are completely disrupted by z = 0. In the fiducial model with ξ_e = 0.033, about 60% of the original sample is disrupted. (Note that this implies that roughly 5 × 10⁷ M_☉ worth of stars in the Galactic stellar halo could be remnants of the disrupted clusters.) With the factor of ξ_e = 0.02, only ∼30% are disrupted. With the factor of ξ_e = 0.01, almost all clusters survive.

The effect of decreasing δ and δ₀ is to shift the peak of the mass function to a lower mass.

We repeated the grid parameter search for the best metallicity distribution for two alternative prescriptions, one with ξ_e = 0.02, δ = δ₀ = 1/3, and the other with ξ_e = 0.033, δ = δ₀ = 1/9. We found that in both cases our model could produce an observationally consistent metallicity distribution. However, lowering either ξ_e or δ significantly alters the mass function away from the data by allowing too many low-mass clusters to survive. Raising ξ_e and δ may improve the mass function, but steers away from recent constraints on the two parameters (Baumgardt & Makino 2003; Gieles & Baumgardt 2008). Therefore, we ultimately conclude that our fiducial prescription (ξ_e = 0.033, δ = δ₀ = 1/3) works best.

For illustration, we list below the properties of the best models in the two alternative prescriptions. For ξ_e = 0.02, δ = δ₀ = 1/3, we find a peak model that has a metallicity distribution with P_KS,Z = 69% and mass distribution with P_KS,M = 0.02%. The parameters of this model are p₂ = 4.4, p₃ = 1.4, p₄ = 0, and p₅ = 0.99.

For ξ_e = 0.033, δ = δ₀ = 1/9, we find a peak model that has a metallicity distribution with P_KS,Z = 19% and mass distribution with P_KS,M = 12%. The parameters of this model are p₂ = 4, p₃ = 1.4, p₄ = 0, and p₅ = 0.97. Its metallicity distribution is shown in Figure 23. The overabundance of metal-poor clusters is clear.

Zoom In Zoom Out Reset image size

The latter alternative prescription (δ = 1/9) predicts the disruption time to scale with cluster mass in a manner nearly identical to Gieles & Baumgardt (2008), if we take that all model clusters have a mean galactocentric frequency ω = (2.4 × 10⁷ yr)⁻¹. Based on the orbit calculations discussed in the next section, we find a similar median value of ω for the sample of model clusters in the fiducial model. We also find no correlation between ω and the cluster mass, which confirms our assertion in Section 3 that in such a model the disruption time of the average cluster would scale with mass as ν⁻¹_ev ∝ M^2/3.

Ashman et al. (1994) popularized a mixture modeling code KMM for detecting bimodality in astronomical applications. This code has been widely used for globular cluster studies and can be considered a standard method in the field. The KMM algorithm assumes that an input sample is described by a sum of two Gaussian modes and calculates the likelihood of a given data point belonging to either of the two modes. It also calculates the likelihood ratio test (LRT) as an estimate of the improvement in going from one Gaussian to two Gaussian distributions. However, the LRT obeys a standard χ² statistic only when the two modes have the same width (variance), which may not be satisfied by real datasets. Even though the probability of the LRT can be estimated using bootstrap in principle, in practice the use of the KMM code has been limited to common width modes (the so-called homoscedastic case). Brodie & Strader (2006) and Waters et al. (2009) provide further discussion of KMM.

The KMM method belongs to a general class of GMM algorithms. GMM methods maximize the likelihood of the dataset given all the fitted parameters, using the expectation-maximization (EM) algorithm (e.g., Press et al. 2007). A major simplification, which allows one to derive explicit equations for the maximum likelihood (ML) estimate of the parameters, is that each mode is described by a Gaussian distribution.

For simplicity, and as appropriate for the metallicity distribution, we consider a univariate input dataset. However, the algorithm is fully scalable to multivariate distributions. The likelihood function of a univariate sample x_n is

$$\begin{equation} {\cal L}_K = \prod _n \left(\sum _{k=1}^K \, p_k \, N(x_n | \mu _k,\sigma _k) \right), \end{equation} \tag{ A1 }$$

where

$$\begin{equation} N(x | \mu,\sigma) = {1 \over (2\pi \sigma ^2)^{1/2}} \exp {\left[-{(x-\mu)^2 \over 2\sigma ^2}\right]} \end{equation} \tag{ A2 }$$

is the Gaussian density. The modal fractions are normalized as ∑_kp_k = 1. A unimodal distribution (K = 1) has two independent parameters (μ and σ), whereas a bimodal distribution has five parameters (p₁, μ₁, σ₁, μ₂, σ₂), since p₂ = 1 − p₁.

The power of the GMM method lies in its ability to determine the ML values of the parameters (p_k, μ_k, σ_k). The disadvantage is that the method will always split the dataset into the specified number of modes, K. In order to detect bimodality, it is extremely important to be able to judge whether the bimodal fit is an improvement over the unimodal fit. For this purpose the KMM code uses the LRT test, which appears to be an approximation derived by Wolfe (1971) for the homoscedastic case (σ₁ = σ₂). Define the ratio of the MLs as $\lambda \equiv {\cal L}_{1,\rm max}/{\cal L}_{2,\rm max}$. According to numerical Monte Carlo studies of Wolfe (1971), the statistic −2ln λ approximately obeys the χ² distribution with a number of degrees of freedom equal to "twice the difference between the number of parameters of the two models under comparison, not including the mixing proportions" (McLachlan 1987). This is the χ²₂ distribution in our case. However, the statistic does not apply in the heteroscedastic case σ₁ ≠ σ₂ (it would have been χ²₄). Note that this unusual number of degrees of freedom was found through an empirical approximation. Unfortunately, no exact estimation exists for the goodness of modal split.

Several variations of the method have been suggested in the literature. McLachlan (1987) proposed a parametric bootstrap to test for the number of components. In this method, a test sample is drawn randomly from a unimodal Gaussian distribution with the parameters {μ, σ} best fitting the input sample. The number of objects in the test sample is taken to be the same as in the input sample. The bimodal split is calculated for this test sample using the EM algorithm and the likelihood ratio λ_boot is saved. Repeating the bootstrap a large number of times, we obtain the probability of randomly drawing the ratio as large as that observed in the input sample, λ_obs. If the probability is below a few percent, we reject the null hypothesis that the input sample belongs to a unimodal Gaussian.

However, the parametric bootstrap is not a perfect solution. In the limit of a large number of objects in the input sample, the likelihood function is very sensitive to outliers far from the center of the distribution. Simple measurement errors in the wings of the Gaussian function may cause a unimodal distribution to be rejected, even if it is correct. In other words, GMM is more a test of Gaussianity than of unimodality (see Muthén 2003; Bauer 2007, for more discussion).

Lo et al. (2001) proposed a modified LRT method to test for the true number of components of a Gaussian mixture. The modified statistic must be evaluated numerically, but it still does not address the problem with Gaussian wings. Subsequently, Lo (2008) suggested to use the standard LRT with the parametric bootstrap to test for a heteroscedastic split, and also suggested restricting the ratio of the standard deviations of the two modes to be not less than 0.25 to avoid numerical artifacts. Such a method was recently implemented in globular cluster studies by Waters et al. (2009).

The sensitivity of LRT to the assumption of Gaussian distribution calls for additional, independent tests of bimodality. A useful and intuitive statistic is the separation of the means relative to their widths:

$$\begin{equation} D \equiv {|\mu _1 - \mu _2| \over \big[ (\sigma _1^2 + \sigma _2^2)/2 \big]^{1/2}}. \end{equation} \tag{ A3 }$$

We use the factor $\sqrt{2}$ for consistency with the definition in Ashman et al. (1994), who noted that D > 2 is required for a clean separation between the modes. If the GMM method detects two modes but they are not separated enough (D < 2), then such a split is not meaningful. The power of GMM in this case is counterproductive. A histogram of such a distribution would show no more than two little bumps, which would not be recognizable as distinct populations.

Another simple statistic is the kurtosis of the input distribution. A positive kurtosis corresponds to a sharply peaked distribution, such as the Eiffel Tower. A negative kurtosis corresponds to a flattened distribution, such as a top hat. The sum of two populations, not necessarily Gaussians, is broader than one population and therefore has a significantly negative kurtosis. However, kurt < 0 is a necessary but not sufficient condition of bimodality. A broad unimodal distribution, such as an actual top hat, also has negative kurtosis. Therefore, kurt < 0 is only useful as an additional check to support the results of LRT and the D-value.

In order to provide a more robust measure of the modal split, we have revised and implemented the GMM algorithm independently of the KMM code. We begin with a single run of the EM algorithm to calculate the means and standard deviations assuming a heteroscedastic bimodal distribution. Then we repeat the estimation assuming a unimodal Gaussian case. We take the ratio of the likelihoods λ, the separation D, and the kurtosis as the three statistics of interest. We then estimate the error distribution for the modal parameters using a non-parametric bootstrap (drawing from the input sample with repetitions) of 100 realizations. We also run the parametric bootstrap to assess the confidence level at which a unimodal distribution can be rejected based on each of the three statistics. Its practical application for an input sample of more than 100 objects is limited to about 1000 bootstrap realizations, limiting the confidence level to ∼10⁻³. A sufficiently low probability of each statistic means a unimodal distribution can be rejected in favor of a bimodal distribution. The code also calculates the probability of each data point belonging to either mode.

A sum of two Gaussians with the same variance can sometimes be preferred to the case of different variances because of the one fewer degree of freedom. The choice between the homoscedastic and heteroscedastic cases can be similarly made using LRT, but we feel that it is less important than choosing between a bimodal and unimodal distribution. For comparison with the KMM code, we calculate the homoscedastic split and its approximate probability using the χ²₂ statistic. We also calculate an alternative split into two Gaussian modes with the same mean but different variance. This is a more extended distribution than a single Gaussian, which may be a better fit for a unimodal but non-Gaussian sample.

The steps of our algorithm are summarized below:

1. Calculate {μ, σ} and {p₁, μ₁, σ₁, μ₂, σ₂} using a single EM run.

2. Form three statistics: λ, D, and kurt.

3. Run non-parametric bootstrap to estimate the errors: Δμ_k, Δσ_k, Δp₁, and ΔD.

4. Run parametric bootstrap to estimate the probability of a unimodal distribution, according to λ, D, and kurt.

As a first test of the algorithm, we verified that it reproduces exactly the test output of the KMM code, given the test input provided with the code. We also made two random realizations of a unimodal Gaussian distribution, N(0, 1), with 150 objects and 1500 objects, respectively. The smaller sample has the mean and standard deviation of μ = −0.088 and σ = 0.987, within the intended target given the sample size. Indeed, a non-parametric bootstrap gives Δμ = 0.084 and Δσ = 0.063. The kurtosis of the input sample is kurt = 0.104. A heteroscedastic split gives two peaks, by construction, with μ₁ = −0.483, σ₁ = 1.206 and μ₂ = −0.032, σ₂ = 0.938. However, the split is not statistically significant. The likelihood is improved only by −2ln λ = 0.26 relative to the unimodal case, which gives the probability better than 99% that the input sample is unimodal. The parametric bootstrap gives a similar probability of 96%. The separation of the peaks also leads to the same conclusion: D = 0.42 ± 1.46. The parametric bootstrap probability of drawing such value of D randomly from a unimodal distribution is 87%. The probability of drawing the measured kurtosis is 74%. Thus, all three statistics correctly show that the input distribution is not bimodal. The larger test sample has smaller parameter errors, as expected, but similar significance levels from the parametric bootstrap.

We then apply the GMM algorithm to the sample of observed metallicities of the Galactic globular clusters. A unimodal fit gives μ = −1.298 ± 0.049 and σ = 0.562 ± 0.028, where the errors are calculated with the non-parametric bootstrap. A heteroscedastic split gives μ₁ = −1.608 ± 0.064, σ₁ = 0.317 ± 0.051 and μ₂ = −0.583 ± 0.074, σ₂ = 0.281 ± 0.075. Of the 148 clusters, 103 (or 70%) are in the metal-poor group and 45 (or 30%) are in the metal-rich group. A homoscedastic split gives μ₁ = −1.620 ± 0.037, μ₂ = −0.608 ± 0.055, and σ₁ = σ₂ = 0.303 ± 0.026. In this case, there are 101 metal-poor clusters and 47 metal-rich clusters. In either case, the likelihood improvement in the 1000 parametric bootstrap realizations is never as high as observed, −2ln λ = 27.5. That is, a unimodal distribution is rejected at a confidence level better than 0.1%. The separation of the peaks is also very clear, D = 3.42 ± 0.47. The observed cluster distribution is indeed bimodal.

When we apply the GMM algorithm to our model sample, we also find strong bimodality. In order to keep the sample size similar to the data, we test separately each of the 11 random model realizations of each host halo. The resulting parameters form a set from which we calculate the average parameters. To estimate both the statistical and systematic errors of each parameter, we sum in quadrature its mean standard deviation in the set and the scatter among the realizations. A heteroscedastic split gives μ₁ = −1.543 ± 0.045, σ₁ = 0.318 ± 0.043 and μ₂ = −0.576 ± 0.077, σ₂ = 0.237 ± 0.026. The breakdown is 66% (±4%) of the clusters in the metal-poor group and 34% in the metal-rich group. The average separation of the peaks is very significant, D = 3.44 ± 0.29. None of the 1000 parametric bootstrap sets of a unimodal distribution returned any of the three statistics (λ, D, kurt) as high as those found in the model. Again, a unimodal distribution is rejected at a confidence level better than 0.1%. Thus, our fiducial metallicity distribution is bimodal as well.

A methodology different from the LRT in the Bayesian approach is the odds ratio, or the Bayes factor described in Liddle (2009). The odds ratio is the ratio of the integrals of the likelihood function ${\cal L}$ over all possible ranges of the model parameters with the corresponding normalized distribution functions of the parameters. In the case of one-Gaussian versus two-Gaussians, it is $\int {\cal L}_2 P(p_1) P(\mu _1) P(\sigma _1) P(\mu _2) P(\sigma _2) dp_1 d\mu _1 d\sigma _1 d\mu _2 d\sigma _2 / \int {\cal L}_1 P(\mu) P(\sigma) d\mu d\sigma$. The odds ratio is dimensionless and is greater than 1 if the modal split indeed significantly improves the likelihood ${\cal L}_2$ over ${\cal L}_1$. On the other hand, if the improvement in the likelihood is small, it can be washed out by integration over the additional parameters and the odds ratio becomes less than 1. We have experimented with using the odds ratio as an objective criterion for the goodness of modal split but found that its application depends sensitively on the adopted range of the parameters and their (unknown) distribution functions. Therefore, we have not included it in our GMM code.

A completely independent test of unimodality was proposed by Hartigan & Hartigan (1985). It was first used for globular cluster studies by Gebhardt & Kissler-Patig (1999). The dip test is based on the cumulative distribution of the input sample. The dip statistic is the maximum distance between the cumulative input distribution and the best-fitting unimodal distribution. In some sense, this test is similar to the K-S test but the dip test searches specifically for a flat step in the cumulative distribution function, which corresponds to a "dip" in the histogram representation. The probability of rejecting a unimodal distribution is calculated empirically and tabulated as a function of sample size. We obtained an updated table of the probabilities calculated recently by M. Maechler (http://www.cran.r-project.org/web/packages/diptest).

We have added a driver routine to the original Fortran code of Hartigan & Hartigan (1985). Our code interpolates the probability table for any input sample size up to 5000 objects. Looking just at the significance levels, the dip test appears less powerful than GMM. The dip probability of the observed Galactic sample being bimodal is 90%, whereas the LRT probability is 99.998% and the parametric bootstrap probability is 99.9%. However, the dip test has the benefit of being insensitive to the assumption of Gaussianity and is therefore a true test of modality. It is also much faster to run than the GMM code. We use the dip probability for assessing bimodality of our model realizations in Section 4.4.

Note that we searched for a statistical method that evaluates bimodality using full input information, without having to bin the data. Another method, RMIX, based on a histogram of the input sample, was recently suggested by Wehner et al. (2008).

Source codes of the GMM and dip tests are available from one of the authors (O.G.) upon request.