© ESO, 2012

1. Introduction

Globular cluster (GC) systems exhibit a bimodal optical colour distribution in the majority of luminous and intermediate luminosity early-type galaxies (Ashman & Zepf 1993; Elson & Santiago 1996; Peng et al. 2006). This is widely interpreted as being due to the presence of two old sub-populations ( ≳ 10 Gyr) that differ in metallicity (Brodie & Strader 2006, and references therein). Metallicity bimodality is obvious for the Milky Way (Zinn 1985; Bica et al. 2006) but seems less evident for our spiral neighbour M31 (e.g. Caldwell et al. 2011).

The presence of two peaks in the GC colour distributions has been often taken as an argument in favour of the existence of two major epochs/mechanisms of star formation in the host galaxies. There has been much debate in the literature over the past decades as to what is responsible for this bimodality (see Brodie & Strader 2006). Among the scenarios for GC formation that account for it, one can list Ashman & Zepf (1992), Forbes et al. (1997), Côté et al. (1998), Beasley et al. (2002), Strader et al. (2005) and Rhode et al. (2005). All of the above assumed different formation channels for blue (metal-poor) and red (metal-rich) clusters. However, the properties of these cluster populations are generally not too distinct. For instance, Peng et al. (2008) studying the specific frequencies of 100 early-type galaxies from the Virgo Cluster Survey find basically the same trends when separating the GC systems in red and blue. Peng et al. (2006) analysing the GC systems colour distributions of the same data set find that there is great similarity between the GC peak-metallicity galaxy-mass relations for the two populations. This implies that the conditions of GC formation for metal-poor and metal-rich GCs could not have been too different.

The complexity of the formation histories of early-type galaxies, within the hierarchical merging framework (e.g. Renzini 2006), are not expected to naturally produce the near universality of bimodal GC metallicity distributions. This is actually the case for models of GC formation built upon this hierarchical framework (e.g. Beasley et al. 2002) where bimodality occurs only after introducing a mechanism that artificially truncates the formation of metal-poor GCs. Nonetheless, more recently, Muratov & Gnedin (2010) introduce a scenario in which the formation of the two metallicity sub-populations of GCs may be a natural outcome of the hierarchical theory of galaxy formation in some, although not the entire range of model realisations. Muratov & Gnedin (2010) prescribe the formation of GCs semi-analytically using assembly histories from cosmological simulations combined with observed scaling relations for the amount of cold gas available for star formation.

The interpretation that colour distributions translate directly into metallicity distributions was challenged by Yoon et al. (2006). They demonstrate that non-linear colour-metallicity relations caused by the horizontal branch (HB) morphology may transform a unimodal metallicity distribution into a bimodal optical colour distribution. This issue has been investigated in more detail by Cantiello & Blakeslee (2007), who conclude that combinations of optical and near-infrared (NIR) colours are much less sensitive to this effect. Richtler (2006) had already shown that a flat metallicity distribution can result in a bimodal colour distribution using the Washington photometric system (acknowledging discussions with Boris Dirsch).

A few observational attempts to address whether optical colour bimodality is really representative of metallicity bimodality exist in the literature, providing sometimes ambiguous or conflicting results. For example, Strader et al. (2007) find clear evidence for two metallicity subpopulations in the spectroscopic sample of 47 NGC 4472 GCs from Cohen et al. (2003). For the same set of data, though, the analysis of Cohen et al. do not strongly favour bimodality. Kundu & Zepf (2007) find an optical/NIR (I − H) bimodal distribution for NGC 4486 in a small sample of GCs from one NICMOS/HST pointing combined with HST/WFPC2 data. Moreover, Spitler et al. (2008) presented an optical/mid-infrared analysis using Spitzer Space Telescope for NGC 4594 and NGC 5128. Both galaxies present a clear optical bimodality. While NGC 5128 presents an obvious R −  [3.6]  bimodal distribution, more compatible with the multimodal peaks that the spectroscopic work of Beasley et al. (2008) find (see also Woodley et al. 2010), NGC 4594 has a less clear bimodal distribution in this baseline. In contrast, Alves-Brito et al. (2011) present spectra for over 200 GCs in the Sa NGC 4594 and find a clear bimodal distribution. More recently, Foster et al. (2010, 2011) find very similar Calcium triplet values for red and blue GCs in NGC 1407 and NGC 4494 despite their large colour difference. Since Calcium triplet is a metallicity indicator, similar values for red and blue GCs indicate similar metallicities. One possible explanation given by the authors is the non-linear conversion between colour and metallicity.

Blakeslee et al. (2010) simulate GC populations with both a mass-metallicity relation (blue-tilt) and a non-linear colour-metallicity relation and find bimodal colour distributions with a blue-tilt even though the metallicity distribution appears unimodal.

Finally, the very recent works by (Yoon et al. 2011a,b) show that colours such as (u − g) and (u − z) have significantly less inflected colour metallicity relations than (g − z). They also show that the metallicity distributions obtained from inflected colour metallicity relations are strongly-peaked, unimodal and with a broad metal-poor tail, similar to that of the resolved field stars in nearby elliptical galaxies and those produced by chemical evolution models of galaxies (e.g. Bird et al. 2010).

The purpose of this work is to investigate the nature of optical/NIR colour distributions of different GC systems in several early-type galaxies. We aim to shed light on the true nature of the metallicity distributions of GC systems.

2. Observations and data

The observations and data reduction techniques of the data used in this study are described thoroughly in Chies-Santos et al. (2011a). Here we briefly summarise the procedures applied. A sample of 14 early-type galaxies was imaged in the Ks-band (from now on referred to as only K) with LIRIS at the WHT and combined with archival ACS/HST F475W (~g) and F814LP (~z) images. The galaxies have M_B < −19 and (m − M)    <    32. ACS images were reduced with MULTIDRIZZLE (Koekemoer et al. 2002) and LIRIS images with LIRISDR, in addition to standard IRAF routines.

The GCs were detected and had sizes (R_eff) measured in the ACS images with DAOFIND and ISHAPE (Larsen 1999), respectively. Aperture photometry was performed in the g, z and K bands with PHOT (Stetson 1987). After the following criteria were applied to automatically detected sources, g < 23, 0.5 < (g − z) < 2.0, 1 < R_eff(pc) < 15 a careful visual inspection was performed where obvious non-cluster objects were removed. Also, sources that were too close together in the ACS images and that appeared as one bigger source in the LIRIS images, due to its lower resolution compared to ACS, were flagged. Finally, we caution the reader that the data of NGC 4382 and NGC 4473 were taken at non-photometric conditions.

3. Integrated colours and horizontal branch morphology

In this section we study the effect of the HB in the integrated colours of GCs. The colour of the HB varies abruptly between [Fe/H]  =  −0.6 and  − 0.9. At this [Fe/H] range the HB departs from the red-clump position (Lee et al. 1994). This is argued as the cause of the non-linear behaviour of the colour-metallicity relation (Yoon et al. 2006). Peng et al. (2006) with data from the Milky Way, NGC 4486 and NGC 4472 showed that the empirical transformation from metallicity to the optical (g − z) colour is clearly non-linear. Dirsch et al. (2003) had already called attention to the fact that the colour-metallicity relation is non-linear and that starting from the bluest colours equidistant colour intervals are projected onto progressively larger metallicity intervals.

In Fig. 1 (g − z), (g − K) and (z − K) colour-metallicity relations are shown. For those relations we used 14 Gyr SPoT¹ (Raimondo et al. 2005; Biscardi et al. 2008) and YEPS² models. The models for (g − z) show a clear departure from linearity around [Fe/H]  ~  −0.5; the wiggly feature noticed by Yoon et al. (2006). This wiggle that is also visible in other commonly used colours such as (V − I) and (B − I) (Cantiello & Blakeslee 2007) is generated by the transition from blue to red HBs. This feature is possibly also recognizable in (g − K), although very mildly. Such feature is not visible when considering (z − K). This last colour does however show a non-linearity at intermediate metallicities more pronounced for the SPoT models with respect to the YEPS. This is probably due to an interpolation effect as the point responsible for this, [Fe/H]  ~  −0.4, is right at the transition between blue and red HBs. Nevertheless the sign of the second derivative of the (z − K) − [Fe/H] relation at the non-linearity is opposite from the one for (g − z) −  [Fe/H]. This indicates the different nature of the non-linearities in these 2 colour spaces. If one would remove the [Fe/H]  ~  −0.4 point, the (z − K) model would bear linearity. The SPoT and YEPS models agree well for the (g − z) and (g − K) colour–metallicity relations. For the (z − K) colour–metallicity relation there is more disagreement between the models, with the YEPS models being systematically bluer in this colour than SPOT for a given metallicity. This also happens for (g − K) but to a smaller extent.

Fig. 1 14   Gyr SPoT and YEPS (g − z), (g − K) and (z − K) colour–metallicity relations according to the legend in the middle panel. The black short-dashed-line in the first panel corresponds to the empirical relation to the Peng et al. (2006) data for Milky Way and Virgo Cluster GCs (Blakeslee et al. 2010).

Fig. 2 (z − K) vs. (g − K) for GCs of NGC 4486 and NGC 4649, indicated by different symbols and with error bars in the bottom panel. Note the wavy feature the data presents around (z − K) ~ 2 and (g − K) ~ 3.2. Upper panel: a running median on the data points; middle panel: 14 Gyr SpOT-Teramo and YEPS SSPs with a realistic treatment of the HB morphology (Raimondo et al. 2005; Yoon et al. 2006; Yoon et al. 2011a,b); Bottom panel: Maraston (2005) with blue and red HB at the lowest metallicities and Padova08 14 Gyr SSPs. See Fig. 16 of Chies-Santos et al. (2011b) for a (g − K) vs. (g − z) plot.

In Fig. 2 the (g − K) − (z − K) diagram for the joint NGC 4486 and NGC 4649 GC sample is shown with different simple stellar population (SSP) models over plotted in different panels. Note the wavy feature the data presents around (g − K) ~ 3.2 and (z − K) ~ 2. A running median on the data points is over plotted in the top panel for better visualisation of the wiggly behaviour. The SPoT 14 Gyr SSP (middle panel) is over plotted where this non-linear feature is also present. SPoT models match the observed colour magnitude diagrams of star clusters fairly well over a wide range of ages and chemical compositions (Brocato et al. 2000; Raimondo et al. 2005). In fact, SPoT models are optimised to simulate the HB spread observed in galactic GCs and this may be one point of uncertainty. Sohn et al. (2006) has shown evidence that the most massive NGC 4486 GCs may have different He abundances and thus different HB morphology than galactic GCs. The YEPS models are also shown in Fig. 2. As SPoT, the YEPS models consider the systematic variation in the mean colours of HB stars as a function of metallicity. Even though the SPoT and YEPS tracks in Fig. 2 do not fit the data well in the redder part of the diagram, they present a wavy feature approximately at the same location as the data: the location pinpointing the transition from blue to red HB morphology. Other SSP models, in which the detailed modelling of the metallicity dependent HB morphology of the SPoT/YEPS models is not included, do not show this behaviour (bottom panels for Padova08³ and for Maraston 2005). The SPoT/YEPS models for the other 2-colour combinations also present this non-linear behavior, although it is not possible to see it in the data using (g − K) vs. (g − z) and (g − z) vs. (z − K) as clearly as it can be seen using (g − K) vs. (z − K). We suspect that this could be an effect of scatter as the K_err is dominant. The errors in (g − K) and (z − K) are not actually independent, but both dominated by the K_err. This means that the uncertainties cause the data to scatter more or less diagonally in the plot. Moreover this diagram has a more similar scale than with the other combinations, where only one of the colours contains K. We refer the reader to Fig. 16 of Chies-Santos et al. (2011b) for a plot of (g − K) vs. (g − z).

Fig. 3 (g − z), (g − K) and (z − K) colour distributions for the different GC systems. Over plotted are the bimodal (blue and red peaks in the corresponding colours) distributions returned by KMM and the sum of the blue and red peaks (in green).

By excluding the [Fe/H]  ~  −0.4 point for the SPoT models, the (z − K)–[Fe/H] would be more linear. However, a similar behaviour is found also for the YEPS models at approximately this metallicity. The (g − K) colour presents a near linear relation to metallicity and also has a broader baseline than (z − K).

4. Colour distributions and bimodality tests

In this section we investigate the behaviour of the (g − z), (g − K) and (z − K) colour distributions. They are shown for the different GC systems in Fig. 3. Note that while a great number of systems look bimodal in (g − z) they appear less and less bimodal as one moves on to the colours that should in principle sample less the HB ((g − K) and (z − K)). The bimodality feature is generally the least apparent in the (z − K) histograms. To quantify whether the distributions are better described by bimodal as opposed to unimodal models, the KMM (Ashman et al. 1994) algorithm was applied to the data. The bimodal (blue and red peaks and their sum) distribution estimates returned by the code are shown in Fig. 3. The outputs of KMM for the different colours and galaxies are listed in Table 1. This test utilises the likelihood ratio test (LRT) to estimate the probability that the distribution of a number of data values is better modelled as a sum of two Gaussian distributions than a single Gaussian. However, it obeys a χ² statistics only when the two modes have the same variance (homoscedastic case). P(KMM) indicates the probability of rejecting a unimodal distribution in favour of a bimodal one, with a low probability indicating that a bimodal distribution is preferred.

The GMM (Muratov & Gnedin 2010) code was also applied to the data. GMM is a robust generalisation of the KMM test. Besides modelling a bimodal distribution with different mode variances (heteroscedastic case), it uses three different statistics: LRT(χ²), distance from the mean peaks (DD) and kurtosis (kurt). Moreover, it makes use of bootstrapping to derive uncertainties for the estimators of these different statistics. Finally, GMM outputs the probability that the three statistics give for rejecting a unimodal fit in favour of a bimodal one: P(χ²), P(DD) and P(kurt). From now on, we refer to P(χ²) as P(GMM). The DD statistic is a measure of how meaningful the split between the two peaks is. The kurt statistic is a measure of peakedness of the distribution. A positive kurtosis corresponds to a sharply peaked distribution whereas a negative kurtosis corresponds to a flattened distribution. A negative kurtosis is a necessary but not sufficient condition for bimodality. The DD and kurt are additional checks to support the results of the LRT since the likelihood function is very sensitive to outliers far from the centre of the distribution and may reject a genuine unimodal distribution.

Before running KMM and GMM objects out of the ranges 1 < (g − K) < 5 and 0.5 < (z − K) < 4 were removed from the distributions. These are objects far from the centre of the distributions and are very likely to affect the results of the bimodality tests by e.g. creating artificial peaks which will be interpreted as separate populations by these statistical tests. For NGC 4486, there is no object in this range, while for NGC 4649 three objects drop off the distribution with these restrictions. The outputs of GMM for the different colours and galaxies are listed in Table 2. For the galaxy NGC 4570, due to the low number of clusters, the test did not converge, and these estimates are not shown.

In Fig. 4, histograms of the different probabilities returned by KMM and GMM are shown. The (g − z) probabilities of these tests are concentrated towards 0 rather than 1, indicating the probable rejection of the unimodal distribution. In fact,  ~80% of the (g − z) distributions of P(KMM) and  ~70% of P(GMM) are consistent with  ≤ 0.05. The P(DD) and P(kurt) are  ≤ 0.1 for  ~40% and  ~46% of the systems, respectively. Therefore a bimodal distribution is favoured for approximately half of the cases when the optical colour is used. In contrast, the (g − K) probability values are spread between 0 and 1. Only  ~15% of the cases have P(KMM) and P(GMM) ≤ 0.05. The P(kurt) is  ≤ 0.1 for only  ~8% of the systems, while none have P(DD) ≤ 0.1. Like (g − K), the (z − K) probability values are also spread between 0 and 1. Compared to the former colour, a slightly larger number of systems have P(KMM) and P(GMM) ≤ 0.05 (~21% and  ~38%, respectively). The P(DD) is  ≤0.1 for only  ~8% of the systems, while none have P(kurt) ≤ 0.1. Note however, that generally for (z − K) and for a few cases of (g − K) these tests assign very few objects for one peak compared to the other, while in (g − z) the number assigned to the different peaks is much more equal. In summary, the higher probability values assigned by KMM and GMM for (g − K) and (z − K) attest that bimodality does get less evident in (g − K) and (z − K) compared to (g − z).

5. The cluster rich galaxies and simulations of the colour distributions

The (z − K) and (g − K) colours, where bimodality appears generally less pronounced, are also the colours with the higher photometric uncertainties (Chies-Santos et al. 2011a). In the present section we investigate whether photometric errors might be the reason bimodality is less apparent in such colours.

First we take the three richest GC systems (NGC 4486, NGC 4649 and NGC 4552) to see how the distributions look when samples with restricted photometric errors are considered. The colour distributions for these cluster rich galaxies are shown in Fig. 5 for different ranges of photometric uncertainties: K_err ≤ 0.5 (whole sample), K_err ≤ 0.1 and K_err ≤ 0.05. The outputs of KMM and GMM are shown in Tables 3 and 4 respectively. When applying these restrictions, a greater fraction of blue objects drops off the sample, as compared to red objects. This is due to the colour dependant scatter, as discussed in Chies-Santos et al. (2011b). For the (g − z) distribution, all probability values are significantly low, indicating bimodality both for samples with K_err ≤ 0.1 and K_err ≤ 0.05. For the other colours the situation changes and distinguishing a unimodal distribution from a bimodal one through KMM and GMM becomes very difficult. It can happen that a probability value goes from  ~1 to  ~0 with the restricted sample, e.g. P(GMM) of (g − K) for NGC 4486. Also, it happens that most objects are assigned to one of the peaks in the (g − K) and (z − K) KMM estimates of the K_err ≤ 0.1 and K_err ≤ 0.5 samples (e.g. for NGC 4649). This exercise does not lead to a clear conclusion regarding the presence of bimodality in the optical/NIR colour distributions, mainly due to the low number statistics when running the tests in the samples restricted in photometric uncertainties. As generally assumed, if the (g − z) distribution is clearly bimodal as it is for NGC 4486 one would expect to see bimodal (g − K) and (z − K) distributions in the K_err ≤ 0.05 sample. Instead, the optical/NIR distributions look much less bimodal. However, note that bimodality is somewhat significant in (g − K) for NGC 4486 for the K_err ≤ 0.05 sample, with P(KMM) = 0.035 and P(GMM) < 0.001.

5.1. Simulations including photometric scatter

Here, we test whether bimodality could be blurred in (g − K) and (z − K) due to the high photometric uncertainties in K. As shown in Blakeslee et al. (2010) the most subtle details in the colour metallicity relation matter for the final transformed distribution. Therefore, to avoid unnecessary complications we do not use SSP models to transform the colours but assume linear transformations. A bimodal distribution (in this case, (g − z)) should maintain its form if linearly mapped to other colour spaces. We find the relations between (g − z) and (g − K), and between (g − z) and (z − K) by fitting the best straight line by a procedure that gives more weight to the objects with less observational errors in both the abscissa and the ordinate using the GC systems of NGC 4486 and NGC 4649, as in Chies-Santos et al. (2011b). The relations are given by: $\begin{array}{ccc}\mathrm{(}\mathit{g}\mathrm{-}\mathit{z}\mathrm{)}& \mathrm{=}& \mathrm{0.465}\mathrm{\ast }\mathrm{(}\mathit{g}\mathrm{-}\mathit{K}\mathrm{)}\mathrm{-}\mathrm{0.349}\\ \mathrm{(}\mathit{g}\mathrm{-}\mathit{z}\mathrm{)}& \mathrm{=}& \mathrm{0.746}\mathrm{\ast }\mathrm{(}\mathit{z}\mathrm{-}\mathit{K}\mathrm{)}\mathrm{-}\mathrm{0.387.}\end{array}$We linearly transform the (g − z) distribution to (g − K) and (z − K) using the inverted relations (1) and (2). In the top panels of Fig. 6 the results of this linear transformation are shown for NGC 4486. Recalling the result of Chies-Santos et al. (2011a) that the photometric scatter in K, as measured by PHOT, is underestimated by a factor of 2, we add randomly a realistic K dependent scatter to the (g − K) and (z − K) transformed distributions. The contributions of the g and z errors to these colours are negligible if compared to the K errors. In order to add the error associated to each K-band magnitude to the colours we first determine an analytic relation between K vs. K_err. The middle left panel of Fig. 6 shows the observed K and K_err, the best fit relation as well as its one sigma deviation. The analytic relation is of the form K_err = N × 10^a ∗ K, where N and a are constants that are different for different galaxies. A random number drawn from a gaussian with dispersion corresponding to the measured one sigma deviation of the K vs. K_err relation is calculated. As in Chies-Santos et al. (2011a) it was shown that the photometric errors are underestimated by a factor of 2, we multiply this random number by two and add it to the analytic K_err. Following this procedure we associate a K-band error to each GC according to its magnitude.

One outcome of the randomly sampled realistic scatter is shown in the middle right panel of Fig. 6 as a scatter plot of large black symbols. For comparison, the parent distribution of errors is also shown as a scatter plot of small red dots. In the bottom panels the resulting (g − K) and (z − K) distributions after the scatter is added are shown. A realistic scatter is sampled randomly and added to the linearly transformed distribution 100 times to generate 100 unique distributions. In order to have a sufficient number of clusters for a robust analysis, we do this only for the most GC-rich galaxies that show obvious bimodality in (g − z) and a nearly equal number of clusters in the blue and red peaks: NGC 4486, NGC 4649, NGC 4552 and NGC 4621.

KMM and GMM are run on the outputs of these simulations. In Fig. 7, histograms of the probability values returned by these statistical tests are shown. For NGC 4486 and NGC 4649 a linearly transformed bimodal distribution in (g − z) remains bimodal both in (g − K) and (z − K) even when realistic photometric scatter is added. This is attested by the statistically significant values of P(KMM), P(GMM) and P(kurt) ( ≤ 0.05) for the great majority of the cases. Nonetheless, the separation between the peaks is not as significant, P(DD) > 0.05. However, for NGC 4552 and NGC 4621 the high probability values show that photometric scatter could blur bimodality in (g − K) and (z − K) in many of the outcomes. Moreover, for NGC 4552, KMM does not converge for many cases, probably because of the low number of clusters. Therefore, the only galaxies for which we can robustly analyse the optical/NIR colour distributions are NGC 4486 and NGC 4649.

This analysis shows that if unimodality is preferred over bimodality in (g − K) and (z − K) it is not an effect of photometric errors in NGC 4486 and NGC 4649. KMM and GMM tests show that their (g − z) distributions are genuinely bimodal. The (g − K) distribution is unimodal for NGC 4486 while bimodal for NGC 4649. For the reddest of the colours, (z − K), bimodality is not significant for any of the galaxies.

5.2. Simulations accounting for scatter in the colour–colour diagrams

The analysis performed above is intended to reproduce the observed K − K_err parameter space, and showed that the photometric uncertainty does not play a major role in “smoothing” the observed (g − K) and (z − K) colour distributions. However, it does not take into account the scatter in the (g − z) vs. (g − K) and (g − z) vs. (z − K) diagrams, that might be an intrinsic scatter, at least in part.

In order to take into account the effect of the scatter in the colour-colour diagrams we performed two more sets of simulations for NGC 4486 and NGC 4649. For the second set of simulations, as before, we linearly transformed the (g − z) colour to (g − K) and (z − K) using Eqs. (1) and (2). Following this, we calculated the residual for each GC from the transformation relations (1) and (2). For each galaxy we then simulated 100 colour distributions (both for (g − K) and (z − K) by adding errors randomly pulled from the array of residuals with replacement to the transformed (g − K) and (z − K) disitributions.

In the second paper of this series (see Fig. 14, Chies-Santos et al. 2011b) it is shown that the scatter estimated through a broken fit in the colour-colour diagrams (rather than a linear fit) is almost entirely consistent with the photometric scatter. However, there we also note that there is still extra scatter found at bright magnitudes that may be attributed to age spreads or other factors. Following this, we performed a third set of simulations that consisted in using broken fit relations for transforming the (g − z) colour to (g − K) and (z − K). These are given by Eqs. (3), (4), (5) and (6). Equations (3) and (5) are valid for (g − z) < 1.187 and Eqs. (4) and (6) for (g − z) > 1.187, $\begin{array}{ccc}\mathrm{(}\mathit{g}\mathrm{-}\mathit{z}{\mathrm{)}}_{\mathit{b}}\mathrm{=}\mathrm{0.260}\mathrm{\ast }\mathrm{(}\mathit{g}\mathrm{-}\mathit{K}{\mathrm{)}}_{\mathit{b}}\mathrm{+}\mathrm{0.232}& & \\ \mathrm{(}\mathit{g}\mathrm{-}\mathit{z}{\mathrm{)}}_{\mathit{r}}\mathrm{=}\mathrm{0.340}\mathrm{\ast }\mathrm{(}\mathit{g}\mathrm{-}\mathit{K}{\mathrm{)}}_{\mathit{r}}\mathrm{+}\mathrm{0.140}& & \\ \mathrm{(}\mathit{g}\mathrm{-}\mathit{z}{\mathrm{)}}_{\mathit{b}}\mathrm{=}\mathrm{0.292}\mathrm{\ast }\mathrm{(}\mathit{z}\mathrm{-}\mathit{K}{\mathrm{)}}_{\mathit{b}}\mathrm{+}\mathrm{0.431}& & \\ \mathrm{(}\mathit{g}\mathrm{-}\mathit{z}{\mathrm{)}}_{\mathit{r}}\mathrm{=}\mathrm{0.437}\mathrm{\ast }\mathrm{(}\mathit{z}\mathrm{-}\mathit{K}{\mathrm{)}}_{\mathit{r}}\mathrm{+}\mathrm{0.391.}& & \end{array}$Similarly to the first set, we ran KMM and GMM on the output of these set of simulations. Moreover, before running KMM and GMM we applied the same colour cuts (1 < (g − K) < 5 and 0.5 < (z − K) < 4) as for the real data. In Fig. 8 we show the outputs of KMM and GMM of the third set of simulations, similarly to Fig. 7. For the third set of simulations these tests converged for  ≳96% of the cases for both KMM and GMM for both simulated colours. The (z − K) simulated distributions of NGC 4486 converged for  ~88% of the cases. Figure 8 only shows the cases which converged. From the left panels of Fig. 8 it clear that there is a peak at 0.05 for P(KMM), P(GMM) and P(kurt). This indicated that for the (g − K) simulations of NGC 4486 and NGC 4649 bimodality is favoured over unimodaliy. Formally,  ~44% of the simulations are found to be significantly bimodal considering P(KMM) for these GC systems. If one considers P(GMM) these values are 38% for NGC 4486 and 32% for NGC 4649.

Fig. 4 (g − z) (top panel), (g − K) (middle panel) and (z − K) (bottom panel) probability values from KMM and GMM.

Fig. 5 NGC 4486, NGC 4649, and NGC 4552 colour distributions in (g − z), (g − k) and (z − k). The open histogram are the colour distributions for all clusters that make the final sample, according to the criteria outlined is Sect. 2. The hashed histograms show the same distributions when only the clusters with Kerr ≤ 0.1 are left in the sample and the filled histogram when only the clusters with Kerr ≤ 0.05 are present.

However, the interpretation of the right panels of Fig. 8 is slightly more tricky. First, the P(kurt) values are skewed to 1 especially for NGC 4486, which attests the strongly peaked unimodal distributions. This, combined with tails of outliers and/or short baseline, will cause KMM and GMM to assign a high probability of bimodality to many of the simulated distributions. In fact, KMM (GMM) assigns a high probability of bimodality (≲0.05) for  ~52%(27%) and  ~42%(17%) of the simulations respectively. KMM and GMM assign the great majority of the objects to one of the peaks and still give high probability of bimodality for several of the simulated (z − K) distributions. However, the parent (g − z) distributions show a nearly equal division between blue and red clusters.

In summary, for the (g − K) distributions of NGC 4486 and NGC 4649 it is reasonable to say that a bimodal (g − z) distribution when transformed to (g − K) will remain bimodal for a little less than half of the cases and at the same time reproduce the observed scatter in (g − z) vs. (g − K). However, one caveat to keep in mind is that the simulations based on the two colour diagrams are a worst-case scenario that assumes all the additional scatter comes from (g − K). The estimate as to how likely it is for a bimodal distribution to remain bimodal in the optical/NIR colours in the presence of a scatter that reproduces the observed spread in the colour-colour diagrams is uncertain for (z − K). We do not show the results of the second set of simulations as they are very similar to the third dataset.

One might ask how well do the positions of the peaks of the (g − K) and (z − K) distributions agree with the (g − z) ones when “transformed back” to (g − z) using relations (1) to (6). For example, for NGC 4486 and NGC 4649, using the GMM values of Table 2 we find that the differences between the “transformed back” and the real distributions is very small, of the order of  ~0.01 for relation (1), involving (g − K) and (g − z); but can be as high as  ~0.2 for relation (2) involving (z − K) − (g − z). When the broken fit relations (3)–(6) are used the differences are  ≲ 0.06, being always higher for the relations involving (z − K).

Given that (g − K) is a good metallicity proxy, if NGC 4486 has a unimodal (g − K) distribution, this analysis suggests an underlying unimodal [Fe/H] distribution in the absence of no scatter in the colour-colour diagrams. This is contrary to NGC 4649, whose (g − K) colour and consequently its underlying [Fe/H], appear genuinely bimodal. If colour bimodality is a consequence of metallicity bimodality as typically assumed, then one would have expected more prominent optical/NIR bimodal distributions compared to the optical colour distributions, because the former trace metallicity better than the latter. However, bimodality may be blurred by scatter and is detected in only half of the realizations.

Fig. 6 Illustration for NGC 4486 on how the simulations with realistic photometric errors were performed. Top panels: the (g − K) and (z − K) colour distributions transformed linearly from (g − z) using relations (1) and (2) respectively. Middle panels: the data and the modelled photometric scatter (Kerr) as a function of K and one outcome of the randomly sampled photometric scatter (Kerrran) as a function of K. Bottom panels: the (g − K) and (z − K) simulations with Kerrran added.

Fig. 7 Histograms of the probability values returned by KMM and GMM on the output of the simulations for the cluster rich galaxies with obvious (g − z) bimodality and nearly equal number of clusters assigned for the blue and red peaks; for (g − K) (left panels) and (z − K) (right panels).

Fig. 8 Histograms of the probability values returned by KMM and GMM on the output of the simulations for NGC 4486 (top panels) and NGC 4649 (bottom panels) for (g − K) and (z − K).

6. Discussion

This study shows the appearance of the optical/NIR colour distributions of GC systems for over a dozen galaxies for the first time. The bimodality feature, clear in purely optical colours, is found to be less pronounced in redder colours. For the most GC-rich galaxies, NGC 4486 and NGC 4649, the K-band photometric scatter does not appear to be responsible for the almost disappearance of this, once well-established feature of GC systems. In the absence of any scatter in colour-colour diagrams such as (g − z) vs. (g − K) the non linearities caused by the HB morphology might indeed be the cause of the optical colour bimodality as suggested by Yoon et al. (2006) at least in some of the GC systems. This seems to be the case for NGC 4486, but not for NGC 4649. However, before strong claims are made as to whether certain galaxies have or have not underlying bimodal distributions, data of better quality are necessary.

In some galaxies, blue GCs tend to be redder at brighter magnitudes. This causes the observed blue-tilt (e.g. Strader et al. 2006; Mieske et al. 2006; Harris et al. 2006) feature in the colour magnitude diagrams of certain GC systems. One might argue that the blue-tilt combined with the higher photometric scatter present in the optical/NIR colours in comparison to the purely optical colour could be in part responsible for the disappearance of a clear bimodality feature in optical/NIR colours of GC systems such as NGC 4486 (known to contain the blue-tilt). We argue that there is no reason to think so. By looking at the NGC 4486 histograms in Fig. 5 for clusters with K_err < 0.05, containing only the brightest GCs of the sample, one can see that the (g − z) distribution shows a very clear bimodal feature while the optical/NIR distributions (especially (z − K)) do not show it as clearly. If the blue-tilt were to be (in part) responsible for the non-clear bimodal feature in the optical/NIR colours it would have to also be in the optical one. Moreover, both blue-tilt and bimodality have been detected using (g − z) data. Furthermore this feature cannot be the main factor behind the absence of a clear (g − K) and (z − K) bimodality in NGC 4486 compared to NGC 4649. The GC systems of both NGC 4486 and NGC 4649 have a clear blue-tilt (Strader et al. 2006; Mieske et al. 2006).

Cohen et al. (1998) presented a spectroscopic sample of GCs in NGC 4486. For a marginal preference of bimodality to be detected at the 89% significance level, they have to exclude a tail of very metal rich-GCs. Their [Fe/H] distribution is very narrow, in fact much more than that of the Milky Way, as shown in their Fig. 20.

The fact that Kundu & Zepf (2007) found a clear (I − H) bimodal distribution for NGC 4486 contrasts with our result for unimodality in (g − K). However, the (I − H) −  metallicity relation, shown in Fig. 9 for SPoT and YEPS models show a wavy feature that is responsible for projecting equidistant metallicity intervals into larger colour bins. Could this be argued as the cause for the bimodal (I − H) distribution found by Kundu & Zepf (2007)? Note that the (I − H) dip in the colour distribution of NGC 4486 (Fig. 1 of Kundu & Zepf 2007) occurs between 1.45 − 1.65, exactly where the wavy feature predicted by the SPoT models is. As shown in Sect. 3, the SPoT (g − K)-metallicity relation is far more linear in this regime. Moreover, the field of view of NICMOS, used in Kundu & Zepf (2007) is  ~4 times smaller than that of LIRIS and ACS. Also, the number of clusters shown in the (I − H) distribution of Kundu & Zepf (2007) is roughly only 1/4 of the number in this study and much more centrally concentrated, better sampling the red than the blue sub-population. Therefore the sample LIRIS/ACS sample is significantly different from the NICMOS one.

Fig. 9 SPoT and Yonsei 14 Gyr (I − H) – metallicity relation.

Fig. 10 Simulations of unimodal metallicity distributions with different means for 15 000 (upper left panel) and 80 GCs (bottom left panel) transformed to (I − H) according to the 14 Gyr SPoT – (I − H) – metallicity relation. The resulting respective (I − H) distributions are shown in the top and bottom right panels.

To further investigate what could cause the bimodality seen in Kundu & Zepf (2007) we performed some simulations of unimodal metallicity distributions and transformed them to (I − H) using the 14 Gyr SPoT model. This is done in order to show that the (I − H) distribution might arise due to this specific colour-metallicity relation combined with the fact the sample is biased to metal-rich GCs. In the top panel of Fig. 10, examples of these simulations of unimodal metallicity distributions with different means for 15 000 (estimated NGC 4486 total number of GCs, Peng et al. 2006) and 80 (number of GCs in the Kundu & Zepf 2007 sample) GCs are shown. They are different gaussian distributions with the same dispersion   =  0.5 and means   =  0.2, 0.7 and 1.2. The resulting respective colour distributions for these representative examples are shown in the bottom panel of Fig. 10. Note that the most metal-rich unimodal metallicity distribution becomes bimodal in (I − H) both for the 15 000 and for the 80 GCs examples. Note also that the dips occur at the same (I − H) values (~1.5) of the distribution of Kundu & Zepf (2007). While the simulation that is more closely related to the NGC 4486 observed metallicity distribution (Cohen et al. 1998) is the intermediate metallicity one, shown in black, the one shown in red matches better the observed data of Kundu & Zepf (2007). This is probably because the (I − H) observed distribution is biased to more metal-rich values. It is well known that GCs in the centres of galaxies are more metal-rich than in their outskirts (e.g. Liu et al. 2011). One might argue that the (I − H) transformations shown in Fig. 10 are idealised as they do not contain photometric scatter. The uncertainty in (I − H) for the Kundu & Zepf (2007) study is found to be less than 0.1. Following this, a 0.1 gaussian scatter was added to the transformed (I − H) distributions. For the most metal-rich distributions (red histograms), bimodality continues to be seen in (I − H) for the great majority of the cases (with and without the addition of a 0.1 Gaussian scatter).

The near universality of GC colour bimodality, once well consolidated as direct implication for metallicity bimodality appears to go contrary to the evidence presented here. In the absence of any scatter in the colour-colour diagrams our analysis shows that the equivalence of bimodal colour and bimodal metallicity distribution is questionable and there is room for some speculation. Not all luminous galaxies necessarily have bimodal metallicity distributions. Each galaxy is particular and full of peculiarities, and its metallicity distribution should account for that. If in the Milky Way [Fe/H] bimodality is absolutely evident, this shows that in our Galaxy the major star formation episodes occurred perhaps in less complicated ways. It might be that NGC 4649 due to its obvious optical/NIR colour bimodality had a similar history to that of the Milky Way. CD galaxies, such as NGC 4486 and NGC 1399 (see the unimodal (I − H) distribution for NGC 1399 found by Blakeslee et al. 2012), are much more prone to interactions with their neighbours due to their physical position in the potential well of the galaxy cluster. Perhaps the quantity and/or strength of such interactions play an important role in determining the shape of the [Fe/H] distributions observed today. This picture of non-universality of [Fe/H] bimodality fits much better with the hierarchical merging paradigm.

High S/N NIR/optical imaging in an 8–10 m class telescope is of the utter most importance to further investigate this issue. Also, large spectroscopic datasets coupled with simulations for several GC systems would shed more light in the true nature of the metallicity distributions of GC systems in large galaxies.

7. Summary and conclusions

We have analysed the (g − z), (g − K) and (z − K) colour distributions of GC systems in a sample of 14 E/S0 galaxies. The results are summarised below:

• 1.

  The data presents a non-linear feature around(z − K) ~ 2 and (g − K) ~ 3.2, marking the transition from blue to red HB morphology. SPoT and YEPS models, with a realistic treatment of the HB also show a similar feature. According to these models for old ages the colour metallicity relation does not present the prominent wiggle of Yoon et al. (2006) for the optical/NIR colours: (g − K) and (z − K).

• 2.

  While the great majority of GC systems present an obvious bimodal distribution in (g − z), bimodality is clearly less pronounced in the optical/NIR distributions (g − K) and (z − K). The two most cluster rich galaxies show some remarkable differences. While the GC system of NGC 4486 shows no obvious bimodality in (g − K) when all the GC sample is considered, the GC system of NGC 4649 does. However, if restricted to a brighter sub-sample with small K-band errors (<0.05 mag) the (g − K) distribution of NGC 4486 GCs is better described by a bimodal distribution. Simulations of the (g − K) and (z − K) distributions with realistic K-band errors suggest that the K-band errors cannot be the responsible for the blurring of genuinely bimodal optical/NIR colour distributions in NGC 4486 and NGC 4649. However, when taking into account the extra scatter present in colour–colour diagrams such as (g − K) vs. (g − z), we find that bimodality is indeed likely to be undetectable in over half of the cases for the (g − K) distributions of these two GC systems.

• 3.

  The underlying metallicity distribution of the GC system of NGC 4649 appears to be a genuine case for bimodality. However, for NGC 4486 the situation is less clear and bimodality is detected at a statistically significant level only for the brightest sub-sample of the clusters in (g − K). This bimodality becomes less pronounced when including objects with larger errors, or for the (z − K) colour distributions. In the galaxy, centre of the Virgo Cluster, the argument put forward by Yoon et al. (2006) might contribute for the clear optical colour bimodality. Higher S/N NIR imaging is strongly needed to understand whether the (g − K) lack of bimodality for NGC 4486 could be due to scatter other than photometric. Also it would further constrain the nature of the optical/NIR colour distributions (and the underlying metallicity distributions) of the other, less GC-rich galaxies in the sample.

Acknowledgments

We thank Oleg Gnedin for providing us with the GMM algorithm and the anonymous referee for putting a lot of effort in reviewing this manuscript. His suggestions, significantly improved the paper.

References

All Tables

All Figures

	Fig. 1 14   Gyr SPoT and YEPS (g − z), (g − K) and (z − K) colour–metallicity relations according to the legend in the middle panel. The black short-dashed-line in the first panel corresponds to the empirical relation to the Peng et al. (2006) data for Milky Way and Virgo Cluster GCs (Blakeslee et al. 2010).
In the text

Fig. 2 (z − K) vs. (g − K) for GCs of NGC 4486 and NGC 4649, indicated by different symbols and with error bars in the bottom panel. Note the wavy feature the data presents around (z − K) ~ 2 and (g − K) ~ 3.2. Upper panel: a running median on the data points; middle panel: 14 Gyr SpOT-Teramo and YEPS SSPs with a realistic treatment of the HB morphology (Raimondo et al. 2005; Yoon et al. 2006; Yoon et al. 2011a,b); Bottom panel: Maraston (2005) with blue and red HB at the lowest metallicities and Padova08 14 Gyr SSPs. See Fig. 16 of Chies-Santos et al. (2011b) for a (g − K) vs. (g − z) plot.

In the text

	Fig. 3 (g − z), (g − K) and (z − K) colour distributions for the different GC systems. Over plotted are the bimodal (blue and red peaks in the corresponding colours) distributions returned by KMM and the sum of the blue and red peaks (in green).
In the text

	Fig. 4 (g − z) (top panel), (g − K) (middle panel) and (z − K) (bottom panel) probability values from KMM and GMM.
In the text

Fig. 5 NGC 4486, NGC 4649, and NGC 4552 colour distributions in (g − z), (g − k) and (z − k). The open histogram are the colour distributions for all clusters that make the final sample, according to the criteria outlined is Sect. 2. The hashed histograms show the same distributions when only the clusters with Kerr ≤ 0.1 are left in the sample and the filled histogram when only the clusters with Kerr ≤ 0.05 are present.

In the text

Fig. 6 Illustration for NGC 4486 on how the simulations with realistic photometric errors were performed. Top panels: the (g − K) and (z − K) colour distributions transformed linearly from (g − z) using relations (1) and (2) respectively. Middle panels: the data and the modelled photometric scatter (Kerr) as a function of K and one outcome of the randomly sampled photometric scatter (Kerrran) as a function of K. Bottom panels: the (g − K) and (z − K) simulations with Kerrran added.

In the text

	Fig. 7 Histograms of the probability values returned by KMM and GMM on the output of the simulations for the cluster rich galaxies with obvious (g − z) bimodality and nearly equal number of clusters assigned for the blue and red peaks; for (g − K) (left panels) and (z − K) (right panels).
In the text

	Fig. 8 Histograms of the probability values returned by KMM and GMM on the output of the simulations for NGC 4486 (top panels) and NGC 4649 (bottom panels) for (g − K) and (z − K).
In the text

	Fig. 9 SPoT and Yonsei 14 Gyr (I − H) – metallicity relation.
In the text

	Fig. 10 Simulations of unimodal metallicity distributions with different means for 15 000 (upper left panel) and 80 GCs (bottom left panel) transformed to (I − H) according to the 14 Gyr SPoT – (I − H) – metallicity relation. The resulting respective (I − H) distributions are shown in the top and bottom right panels.
In the text