THE GLOBULAR CLUSTER SYSTEM IN M87: A WIDE-FIELD STUDY WITH CFHT/MEGACAM^*

By far the largest collection of globular clusters (GCs) within any one galaxy in the local universe is in M87, the cD giant in the Virgo cluster. Its total GC population is N ≃ 14000 (McLaughlin et al. 1994; Tamura et al. 2006), roughly a third of the GCs population in all the Virgo galaxies combined (see the Virgo Cluster Survey results of Peng et al. 2006). In addition, because both the foreground reddening and field contamination are low, the M87 GCS can be traced outward to extremely large projected galactocentric distance, 100 kpc or more (Tamura et al. 2006). This combination of characteristics makes M87 a unique testbed for systematics of GC systems.

A particular topic of recent interest is the correlation between GC metallicity and luminosity (or mass), a trend which has been found to show up on the blue, metal-poor sequence (Harris et al. 2006; Strader et al. 2006; Mieske et al. 2006; Spitler et al. 2006; Forte et al. 2007; Wehner et al. 2008; Bassino et al. 2008; Harris 2009, and references cited there). In brief, the mean GC metallicity along the blue sequence is seen to increase steadily, though perhaps not linearly, with increasing cluster mass. So far, the most likely physical explanation that has been advanced for the origin of this effect is some form of self-enrichment during the protocluster epoch (Harris et al. 2006; Mieske et al. 2006; Strader & Smith 2008; Bailin & Harris 2009). The essential idea of this interpretation is that any self-enrichment should be more efficient within the more massive, deeper potential wells of the larger clusters.

M87 has turned into something of a flashpoint for discussion of this mass/metallicity relation (MMR). Strader et al. (2006) first claimed that the mean GC color along the entire blue sequence increased roughly linearly with absolute magnitude, which would correspond to a simple power-law scaling of heavy-element abundance with cluster mass, Z ∼ M^p, for a constant mass-to-light ratio. Their measurements consisted of simple fixed-aperture photometry of ∼1000 GCs around M87 from the single-orbit exposures in (g', z') of the galaxy that were taken as part of the Virgo Cluster Survey (Côté et al. 2004). Mieske et al. (2006) used the same material in the full Virgo Cluster Survey analysis in which the cluster photometry was done with point-spread function-convolved (PSF-convolved) King profiles. They found much the same trend, though it was less clear that the MMR was as straightforward as a linear fit in the color-magnitude plane. Forte et al. (2007) used ground-based C, T₁ photometry with standard PSF-fitting photometry, and found a roughly linear correlation of (C − T₁) color versus magnitude along the blue sequence. Waters et al. (2009), on other hand, find no net MMR along the blue sequence (that is, they find a vertical sequence of magnitude versus color) by using size-corrected aperture photometry of about 2000 clusters from deeper HST/ACS images of M87 in (V, I). Most recently, Peng et al. (2009) have analyzed the same ultra-deep HST/ACS data with the PSF-convolved King profile technique described in Jordán et al. (2009) to find that a significant color/magnitude slope exists for the upper half of the GC magnitude range (roughly, for the clusters brighter than the turnover point of their luminosity function). Thus five different measurement procedures from four partially different data sets have yielded results that are largely in mutual agreement, but also show points of disagreement that are still of some importance and need to be settled.

For comparison, the largest GC database used specifically to analyze the MMR is the sample of ∼10, 000 clusters in six giant Brightest Cluster Galaxies (BCGs; Harris 2009). There, the photometry of the individual GCs is obtained from PSF-convolved King-model profiles and is taken in the metallicity-sensitive (B − I) index and with deep exposures.¹ For clusters brighter than the GCLF turnover the data have high S/N, and the total sample is so large that the GCLF can be traced upward to the most luminous known GCs at L ≳ 10⁷L_☉. Extensive evidence is presented that the MMR has a nonlinear form: the blue GC sequence is nearly vertical (i.e., no trend of GC color versus luminosity) for L ≲ 5 × 10⁵L_☉ and then curves toward redder color at higher luminosity.

The quantitative self-enrichment model of Bailin & Harris (2009) reproduces this nonlinear MMR form fairly accurately. The key physical point in the model is that protoclusters less massive than about 10⁶M_☉ do not have deep enough potential wells to hold back much of the first round of SNII-enriched gas, so for these lower mass GCs there should be no MMR on either the blue or red sequence. Conversely, for the higher mass GCs with their deep potential wells, self-enrichment is almost automatically required. A unique prediction of this particular model is that both the red and blue GC sequences should exhibit a MMR at the top end, though it should be less noticeable for the red sequence because it starts from a much higher level of pre-enrichment.

Settling the issue on observational grounds of the detailed form of the MMR (and indeed whether or not it exists in all galaxies) is of some urgency because confusion over the form of the correlation has already arisen in the literature. Yet it is likely to bear quite directly on understanding the formation of massive GCs and their early star formation history (Harris et al. 2006; Mieske et al. 2006; Strader & Smith 2008; Bailin & Harris 2009).

In this paper, I discuss new photometry of the M87 GCs obtained with Canada-France-Hawaii Telescope/Megacam (CFHT/Megacam). The data comprise a large enough GC sample (N ≃ 10⁴) to analyze the MMR to deep limits in addition to correlation with galactocentric distance. The sample size provides the obvious advantage that features such as the MMR, the spatial distribution, or the GCLF can be studied within a single galaxy without having to put together composite samples. The results discussed in the following sections show well-defined blue and red GC sequences extending to high mass, and exhibit a nonlinear MMR quite similar to those in the six-galaxy BCG sample. Though M87 alone can hardly settle the question of the universality of the MMR, these results reinforce the need to investigate this intriguing effect further. In the present paper, I concentrate on the various features of the color and metallicity distributions.

In this discussion, I adopt a distance modulus (m − M)₀ = 31.0 for M87 and the core of the Virgo cluster (Jerjen et al. 2004; Caldwell 2006; Williams et al. 2007) and a foreground absorption and reddening A_I = 0.04, E_{V − I} = 0.03. In the Sloan Digital Sky Survey (SDSS) filter system the corresponding values relevant here are A_i = 0.04, E_{g − i} = 0.04, E_{g − r} = 0.02.

The 36 CCDs making up the CFHT Megacam array provide a 1° square field of view, and the camera is normally used with the u'g'r'i'z' filter set (in what follows I will drop the accents and refer to the magnitudes simply as ugriz). The image scale is 0187 pixel⁻¹, and the entire field makes up a 20 K × 20 K array. The data used in this paper consist of a single set of homogeneous exposures centered on M87 taken in gri during 2004 March and April (Harris, PI), as listed in Table 1. While this material will eventually be surpassed by the much more extensive CFHT/Megacam Next-Generation Virgo Cluster Survey by Ferrarese et al. (see http://www.astrosci.ca/NGVS), much can already be extracted from the earlier set of M87-centered exposures.

Table 1. CFHT/Megacam Images

Filter	N(exp)	Total t (sec)	FWHM
g'	11	3520	074
r'	10	3000	080
i'	11	3300	082

Download table as:  ASCIITypeset image

Individual exposures in each filter were taken with a "small dither" pattern sufficient to splice over the small gaps between the individual CCDs in the array, but the much larger gaps between the top two ranks of CCDs and the bottom two ranks are still present (standard large-dither patterns are available for the camera that would cross these gaps, but these make the astrometric re-registration more difficult). These remaining gaps show up as two ≃350 pixel-wide blank E/W stripes in the array, making up about 3.5% of its total area on the sky. The hour-long total exposure times in each filter were chosen to reach closely similar limiting absolute magnitudes for the globular clusters. The combined images in each filter are those produced through the Megapipe facility of CADC, and were downloaded from the CADC Archive ready for photometry.

The GCs around M87 are seen against a population of field objects including foreground Milky Way stars and background galaxies. Because the Virgo field is at high latitude, the contamination from stars is minimal and the main contaminants are small, faint galaxies. However, much of the field galaxy population can be removed by morphology because the GCs themselves appear star-like on these ground-based images. At the 16 Mpc distance of Virgo, an average GC half-light diameter of 5 pc (Harris 1996) is equivalent to 006, which in turn corresponds to a FWHM of 004 for an average King model concentration log(r_t/r_c) = 1.5 (Larsen 1999). This scale size is thus 20 times smaller than the 08 FWHM seeing on the combined Megacam frames in this data set, making any such objects indistinguishable from stars. Extensive quantitative tests of PSF-convolved globular cluster profiles for distant galaxies such as these show explicitly that the intrinsic profile width falls below the threshold of measurement if its FWHM is less than 10% of the PSF width (e.g., Larsen 1999; Harris 2009). In practice, what this means for the present data is that the GCs can be accurately measured either through conventional fixed aperture or PSF-fitting photometry. Because the crowding levels are minimal, aperture photometry is adopted here (see below), with an aperture radius set at its optimal value similar to the PSF FWHM. Use of aperture photometry also alleviates concerns over minor changes in the PSF shape across the field.

As another basis of comparison, it is worth noting that imaging of the M87 clusters with ground-based image quality of 08 is almost exactly analogous to imaging clusters in the Coma galaxies at d = 100 Mpc with HST image quality of 01. Explicit tests of PSF-convolved GC profiles for the Coma galaxies (Harris et al. 2009) with deep HST data verify that the intrinsic sizes of the GCs are too small to affect their integrated magnitudes from either PSF-fitting or fixed-aperture photometry at levels larger than ≳0.02 mag. Even further, when these second-order magnitude corrections are differenced between two filters (g − i or g − r in this case), the effects of aperture size corrections on color are even smaller third-order ones that can be completely neglected. Similar comments apply to the (C − T₁) ground-based photometry of M87 and NGC 1399 by Forte et al. (2007) in which the seeing was 10 − 16.

The first step in the photometry was to find candidate objects on all three master images and to eliminate the majority of nonstellar objects. The SExtractor code (Bertin & Arnouts 1996) was run to do this, and plots were generated of the SE effective radius and FWHM versus aperture magnitude, as shown in Figure 1. Objects falling outside the generous "stellar" boxes marked in the figures are either very bright, saturated stars—obviously nonstellar extended galaxies—or extremely faint objects, all of which could be eliminated from the detection lists.

Zoom In Zoom Out Reset image size

Next, the SE-culled finding lists defined separately in g, r, i were compared, and any objects found on all three filters and coordinate matched to within ±1.5 pixel were kept for the final photometric measurement run. These matched objects were measured through daophot with an aperture radius of 3.5 pixel (065, equivalent to a linear size of 50 pc at the distance of M87). Crowding in this high-latitude field is negligible, and comparisons of the fixed-aperture magnitudes with PSF-fitted magnitudes showed no differences in the resulting color–magnitude diagrams.

The Megacam-stacked images in the CADC Archive have internally homogeneous astrometry that is already close to the most commonly used systems (such as SDSS, 2MASS, USNO), but to ensure that the data were on a specific astrometric system as closely as possible, the USNO-B1 system was selected. The USNO-B1 catalog was searched for all objects within the Megacam field, and the positions of these objects were matched up with the (α, δ) values returned directly from the Megacam image headers. This matchup yielded 932 stars in common to within 1''. Small zero-point offsets equal to Δα = (0053 ± 0008), Δδ = (0254 ± 0008) in the sense (USNO-CFHT) were found, but no significant higher order terms. Figure 2 shows the residuals after the native Megacam astrometric coordinates were shifted onto the USNO-B1 system. The object-to-object rms scatter is ±023 across the entire field, with no systematic deviations larger than ±005.

Zoom In Zoom Out Reset image size

The Megacam ugriz magnitudes are known to differ slightly from the SDSS scales (see the Megacam webpages for updated calibration information). However, for the purposes of more direct comparison with other GCS data in the literature, it is useful to put the present M87 data onto the SDSS scale. The SDSS DR5 catalog was searched for objects in the M87 field and compared with the Megacam object lists, with the result that matches were found for about 1000 stars brighter than i = 21.0. These were used directly to set the photometric zero points for gri in the M87 field. The mean differences Δm = m(CFHT) − m(SDSS) were found to be Δ(g) = 0.254 ± 0.004, Δ(r) = 0.407 ± 0.008, and Δ(i) = 0.412 ± 0.006, where the CFHT magnitudes are those defined simply using the zero points in the image headers.

A secondary advantage of tying these data closely to a widely used standard system is that the intrinsic colors of globular clusters in the (g − r) and (g − i) indices have not previously been well established. The Virgo Cluster Survey data are in (g − z) through the HST/ACS filters (Peng et al. 2006; Mieske et al. 2006). Previously published studies of GCSs for various galaxies in (g − i) include Forbes et al. (2004), Wehner et al. (2008), and Cockcroft et al. (2009) all taken with Gemini/GMOS, but in all of these the absolute photometric calibrations are preliminary and not likely to be accurate to better than ±0.1 mag (see these papers for detailed discussion). M87 offers an excellent opportunity of setting up a fiducial color scale for GCs, not just because a direct matchup to the SDSS catalog is possible, but also because its foreground reddening is almost negligible and the statistical size of the GC population allows the mean colors of the red and blue sequences to be defined internally to ±0.02 mag or better (see below).

In Figure 3, the final data for 20189 star-like or near-starlike objects are shown, plotted in (i, g − r) and (i, g − i). The detection rate for objects within the color range (g − r) < 1.0 is highly complete for i < 23.5; in the following analysis, the main conclusions rely even more conservatively on the interval i < 23.0 and are unaffected by any incompleteness.

Zoom In Zoom Out Reset image size

Several clearly distinct components show up. The GCs that we are after lie on the two nearly vertical, narrow sequences at (g − i) ≃ 0.8, 1.05, while the populations bracketing these on the blue and red sides are the expected loci of foreground Milky Way halo and disk stars. Any objects falling in the regions

$$\begin{equation*} \eqalign{ & (g-r) < 0.43, \quad (g-r) > 0.95\\ & (g-i) < 0.59, \quad (g-i) > 1.40\\ & i < 18.5 } \end{equation*}$$

(the boundaries of which are marked by the boxes in Figure 3) are thus rejected as field contamination. This final culling leaves 11946 objects that make up our final "best" GC candidate list.

The distribution of the GC candidate colors by spatial location is shown in Figure 4. Here, the projected distance R from the center of M87 is expressed in units of R_eff = 158, the effective radius of the galaxy's spheroid light. A well known though still striking feature of the M87 GC system is its huge radial extent: the clusters extend detectably beyond even R ∼ 150 kpc or 20 R_eff (see the following section). The first hints of its huge GCS halo were already evident in the wide-field star-count work of Harris & Smith (1976) and more clearly outlined with the CCD-based photometry of McLaughlin et al. (1994) and Tamura et al. (2006). For comparison, however, the V magnitude data of McLaughlin et al. (1994) could be used to trace the GCS out only to the radial limits of their single-CCD field at R(max) ≃ 8' ≃ 40 kpc. Tamura et al. (2006) used a combination of four Subaru Suprime-Cam pointings in three colors (BVI) to make up a surveyed region of 136' × 27', about the same total area as the single Megacam field but oriented along the E/W direction. M87 itself is in the SW corner of their westernmost SCam pointing, while another Virgo giant NGC 4552 is in their third field pointing.

Zoom In Zoom Out Reset image size

An equally striking feature of the spatial distribution from Figure 4 is the difference between the blue and red GC subpopulations. The red sequence virtually disappears into the general spread of field contamination for R ≳ 6 − 8R_eff (∼12' or 50 kpc), while the blue sequence continues detectably outward much further. Two other illustrations of this are shown in Figures 5 and 6. As a first approximation, we define the blue GCs as those in the color interval 0.59 ⩽ (g − i) ⩽ 0.93, and the red GCs as those in 0.93 < (g − i) ⩽ 1.40. In Figure 5, the spatial positions of each type are shown. Tamura et al. (2006) and Forte et al. (2007) find exactly the same phenomenon, and these papers provide a more detailed discussion of the spatial distributions of each component. In Figure 6, the CMDs in (i, g − i) are displayed for five radial zones. Whereas the blue and red sequences are almost equally populated for the innermost zone R < 3', the red sequence dies away more rapidly outward, leaving the blue sequence to dominate at large radii.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The M87 GCS is so populous and extended that the level of background contamination (that is, the field contamination still present in the list of starlike objects, assumed constant across the field) is hard to determine precisely. The contamination is completely unimportant at small radii, but at the largest radii of the survey it determines the effective limit beyond which the GCS can be traced. We have no truly remote "control field" adjacent to M87 for measuring σ_b, and ultimately a definitive measurement of the background field level and the outermost radial profile of the M87 GCS will have to await the still wider-field work from the NGVS. For the present, the outer parts of the single Megacam field itself must be used to gauge σ_b. A reasonable estimate can be made by noting that for R ≳ 12R_eff, the areal density of red objects (i ⩽ 23.0 and 0.93 < (g − i) < 1.4) is nearly uniform at σ = 0.36 arcmin⁻², accounting for the edges of the frame and the two blank stripes in the CCD array. The color range enclosing the blue GC sequence is narrower (0.59 ⩽ (g − i) ⩽ 0.93), so if the color distribution of the contaminants is roughly uniform over this range, the level of the blue background will be proportionately smaller, σ_bkgd(blue) ≃ 0.26 arcmin⁻², for the same magnitude range i < 23.0. (See also Harris 2009 for more extensive tests and discussion of the background color distribution in similar high-latitude fields.) Thus, for the present purposes, I adopt a total σ_bkgd ≃ 0.6 arcmin⁻², with a true external uncertainty likely to be near ±0.1 arcmin⁻².

In the final list of GC candidates defined above, there are 8313 objects brighter than i = 23.0. Subtracting off the background, this leaves a total GC population N_GC ≃ 6200 ± 350 brighter than M_I ≃ −8.5, which is about 0.2 mag brighter than the turnover (peak) point of the luminosity function.

The actual radial profile of both the blue and red GC subsystems can be quantified by tracing the projected surface density σ_cl directly as a function of radius. The observed number density of all the candidate GCs is by definition σ = σ_cl + σ_b. The resulting background-subtracted data for i ⩽ 23.0 are listed in Table 2 and shown in Figure 7. The columns in the Table give the inner and outer radii of each annulus, the mean radius, the calculated total number of clusters in the annular zone after background subtraction, and the density σ_cl. The conventional log/log form (log σ_cl versus log R) is shown in the left panel, indicating that both the blue and red components fall off progressively more steeply with increasing radius. A standard Hubble profile, which is a power law with a core radius R₀,

$$\begin{equation} \sigma _{\rm cl}(R) \, = \, \sigma _0 \left(1 + {R \over R_0} \right)^{-a} \end{equation} \tag{ 1 }$$

provides a reasonable fit to either component. Fits of this model to the data (shown in Figure 7) for the blue GCs give (σ₀ = 66 arcmin⁻², R₀ = 20, a = −1.8), while for the red GCs, (σ₀ = 150 arcmin⁻², R₀ = 12, a = −2.1). These are equivalent to simple power laws σ ∼ R⁻ⁿ where the exponent n steepens continuously with radius

$$\begin{equation} n(R) \, = \, a \, { (R/R_0) \over \big(1+{R \over R_0}\big) } \,. \end{equation} \tag{ 2 }$$

Both components also follow a deVaucouleurs-law form (log σ_cl versus R^1/4) fairly well over their entire observed run. In this form, the best-fit equations describing each one are

$$\begin{equation} \log \sigma _{\rm cl}(blue) = (3.154\pm 0.067) - (1.760 \pm 0.045) (R/R_{\rm eff})^{1/4}, \end{equation} \tag{ 3 }$$

$$\begin{equation} \log \sigma _{\rm cl}(red)=(3.874\pm 0.054) - (2.439 \pm 0.037) (R/R_{\rm eff})^{1/4}, \end{equation} \tag{ 4 }$$

where again the σ's are in units of arcmin⁻² and the slopes are normalized to the R_eff of the halo light. (These lines should, however, not be extrapolated inward where the density distribution flattens off more strongly.) For comparison, Tamura et al. (2006) found slopes (−1.47 ± 0.03)(blue) and (−2.44 ± 0.06)(red). As they discuss more extensively, the integrated light of the M87 halo has a R^1/4-law slope of −2.36 ± 0.09, strongly suggesting that the more metal-rich GC population can be directly linked to the similarly metal-rich population of the field halo stars. Much other discussion of similar results from several other galaxies reinforcing this interpretation and connecting with their evolutionary histories can be found in Lee et al. (1998), Lee et al. (2008), Harris & Harris (2002), Bassino et al. (2006), Forte et al. (2007), and Shin & Kawata (2009) among others, and will not be repeated here.

Zoom In Zoom Out Reset image size

Table 2. Radial Distributions

R Range	R	ncl	σcl (arcmin−2)
Blue GCs
100–130	114	68.4	31.57 ± 3.83
1.30–1.69	1.48	84.1	22.94 ± 2.52
1.69–2.20	1.93	128.4	20.75 ± 1.84
2.20–2.86	2.50	135.3	12.93 ± 1.13
2.86–3.71	3.26	206.4	11.67 ± 0.83
3.71–4.83	4.23	256.2	8.57 ± 0.55
4.83–6.27	5.50	293.9	5.82 ± 0.36
6.27–8.16	7.15	388.8	4.56 ± 0.26
8.16–10.60	9.30	428.5	2.97 ± 0.18
10.60–13.79	12.09	423.6	1.74 ± 0.14
13.79–17.92	15.72	426.8	1.04 ± 0.12
17.92–23.30	20.43	472.9	0.68 ± 0.11
23.30–30.29	26.56	382.1	0.33 ± 0.10
Red GCs
100–130	114	81.2	37.47 ± 4.18
1.30–1.69	1.48	102.7	28.03 ± 2.79
1.69–2.20	1.93	135.8	21.93 ± 1.90
2.20–2.86	2.50	143.2	13.69 ± 1.16
2.86–3.71	3.26	167.6	9.48 ± 0.75
3.71–4.83	4.23	170.2	5.70 ± 0.46
4.83–6.27	5.50	186.8	3.70 ± 0.30
6.27–8.16	7.15	189.3	2.22 ± 0.20
8.16–10.60	9.30	187.1	1.30 ± 0.15
10.60–13.79	12.09	166.2	0.68 ± 0.12
13.79–17.92	15.72	145.7	0.35 ± 0.11
17.92–23.30	20.43	94.4	0.14 ± 0.10

Download table as:  ASCIITypeset image

The behavior of the deVaucouleurs-type profile fit itself, near its outer edge, provides a useful consistency check of our adopted backgrounds σ_b(blue,red). If we have chosen the background nearly correctly, then the resulting σ_cl profile can plausibly be claimed (as is the case for most other elliptical galaxies) to follow a nearly straight line in the log σ_cl versus r^1/4 plane all the way to the outer edge of the data. An overestimate of σ_b would cause the background-subtracted profile to drop suddenly and unrealistically to zero at the edge, whereas an underestimate would cause σ_cl to flatten off artificially at a large radius. Neither one is the case (Figure 7b), so the choice of background levels appears to be self-consistent within their stated internal uncertainties.

The σ_cl data indicate that the blue GCs exist in detectable numbers out to at least R = 28' or ≃130 kpc. A final check of the present results comes from comparison with Tamura et al. (2006), whose SCam data extend along a line towards NGC 4552 out to 120' (more than 500 kpc) from M87. They find the blue-sequence GCs to extend detectably to R = 30' (see their Figures 5 and 7), and present a model fit to the combined M87/N4552 counts that suggests the M87 GCs could even extend to 100' and more. Clearly, we have not yet traced the full extent of this system.

A second-order feature of the radial color distributions (Figure 4) is a progressive change in the mean color of both the red and blue GCs with radius: in other words, both components show a metallicity gradient. This feature is clearest for R ≲ 8R_eff or about 60 kpc; beyond that point, the red clusters are almost absent and the blue clusters show no significant change in mean color. For the region inward of 8R_eff, power-law fits yield

$$\begin{equation} \langle g-i \rangle (blue) = 0.813 - 0.025 \log (R/R_{\rm eff}), \end{equation} \tag{ 5 }$$

$$\begin{equation} \langle g-i \rangle (red) = 1.095 - 0.035 \log (R/R_{\rm eff}), \end{equation} \tag{ 6 }$$

where both slopes are uncertain by 10%. The single-stellar-population models of Maraston (2005), which provide integrated color indices for a variety of broadband colors including the SDSS filters, show that Δ(g − i)/Δ[Z/H] = 0.21 ± 0.05 over the metallicity range [Z/H] ≲−0.2 that covers most GCs. The color slopes then convert into metallicity gradients Z ∼ R^−γ where γ = 0.12 ± 0.02 (blue) and 0.17 ± 0.03 (red).

Comparable data for the GCSs in several other large galaxies are now in the literature. Geisler et al. (1996) found γ = 0.15 ± 0.03 (blue), 0.12 ± 0.06 (red) in the other Virgo supergiant NGC 4472. Forte et al. (2001) found γ ∼ 0.2 for the blue GCs in two Fornax ellipticals, NGC 1399 and 1427. Most recently, Harris (2009) derived a mean γ = 0.10 ± 0.02 for both the blue and red GCs in six giant BCGs. All these results are similar, and begin to present a very consistent picture: the mean GC metallicity in these large galaxies depends on spatial location, becoming more enriched the further inward it resides in the halo.

The widespread presence of metallicity gradients, albeit shallow, indicates that both the major epochs of GC formation, metal-poor and metal-rich, had a significant component of dissipative collapse of protogalactic gas and in situ formation which led to higher enrichment levels deeper in the potential wells of these large galaxies. Although gradients in the range γ ∼ 0.1–0.4 for the dominant metal-rich field-star component have been found in many ellipticals within their central bulges (R ≲ R_eff), the globular cluster data indicate that dissipative formation extended much further outward into their halos, to ∼5R_eff and beyond. What is perhaps equally important for understanding galaxy formation is that the metal-poor component, which is all but undetectable in the integrated field-star light, also has an in-built metallicity gradient to the same degree. This result makes it unlikely that the metal-poor stars all formed in pregalactic dwarfs that only later began to merge into bigger galaxies. Instead, the first major metal-poor round of star formation must have taken place while the initial stages of hierarchical merging were already underway (see Harris 2009 for further discussion and references).

We now look more closely at the blue and red GC sequences in the color–magnitude diagrams. The goal is to characterize each of the two sequences through their mean colors as a function of magnitude (the MMR) and their internal color (metallicity) dispersions. To define the MMR (or lack of one) we adopt the same basic approach as used in several previous studies (e.g., Harris et al. 2006; Mieske et al. 2006; Wehner et al. 2008; Bassino et al. 2008; Harris 2009, among others): the CMD is divided into narrow magnitude bins, and within each bin the color distribution of the GCs is fitted by a bimodal Gaussian model. Because each magnitude bin is treated independently, no a priori assumptions are made about the particular form of the MMR, whether linear or otherwise.

For the multimodal fitting, the code RMIX is used here (for extensive descriptions of RMIX applied to this situation, see Wehner et al. 2008 and Harris 2009). The default approach—facilitated by the sheer statistical size of the M87 GC population—is to allow the code to solve freely for all five basic parameters in a Gaussian bimodal distribution: the mean colors of the two modes μ₁, μ₂ (blue, red), their dispersions σ₁, σ₂ (blue, red), and the relative numbers (proportions) p₁/p₂ of objects in each mode. The MMR along either of the sequences is then defined by any trend in the mean μ with magnitude.

The bimodal fits were done for objects within the zone R < 6.5R_eff, i.e., within the part of the M87 halo where both blue and red sequences are definitely present and in roughly equal numbers (Figure 4). We restrict the sample to objects in the color range 0.65 < (g − i) < 1.30, and as a preliminary step, the shallow metallicity gradient described in the previous section is removed by normalizing the color indices to the values that all the clusters would have if they were at R = R_eff. That is, for the blue-sequence GCs (defined as those objects with (g − i) < 0.95), we define (g − i)_norm = (g − i) + 0.0245log(R/R_eff), and for the red-sequence GCs ((g − i)>0.95), we define (g − i)_norm = (g − i) + 0.0345log(R/R_eff).

The results of these bimodal fits in both CMDs are summarized in Tables 2–4. Four separate runs were done with the following conditions, which were chosen to allow additional tests of the reality of any deduced MMR along both of the sequences.

• 1.  
  The data were divided into 0.25 mag bins in i, and within each bin all five free parameters (μ₁, μ₂, σ₁, σ₂, p₁/p₂) defining the bimodal Gaussian fit were freely determined by the fit (i.e., unconstrained). The exceptions to this were at the bright end, where for i < 20.5 half-magnitude bins were used to capture enough objects for a statistically significant fit. In addition, preliminary runs showed that for the four brightest bins (i < 21.0) the fully unconstrained fits tended to be unstable, so for these, the intrinsic widths of the two sequences were constrained to be equal to their average values on the lower bins, σ₁ ≡ 0.060 and σ₂ ≡ 0.10. The results for this run are summarized in Table 3. The columns list the i-magnitude range defining each bin, the number of objects in the bin, the best-fit blue-sequence parameters (μ₁, σ₁) and their uncertainties, the red-sequence parameters (μ₂, σ₂) and their uncertainties, and the best-fit proportion p₁ (the fraction of the total population belonging to the blue sequence). By definition p₂ = 1 − p₁.
• 2.  
  The same conditions as above were used, except that the red sequence was constrained to have a mean color μ₂ ≡ 1.070, the average obtained from the lower luminosity bins in the first (unconstrained) run. The color of the blue sequence was determined freely from the fit. The results for this run are summarized in Table 4.
• 3.  
  The same conditions were used, except that next the blue sequence was pinned at a mean color μ₁ ≡ 0.80, the average obtained from the lower luminosity bins in the first (unconstrained) run. The color of the red sequence was determined freely from the fit. The results for this run are summarized in Table 5.
• 4.  
  The mean colors of both sequences were pinned (μ₁ ≡ 0.80, μ₂ ≡ 1.07) and only the σ's and proportions were solved for.

As an illustration of the overall quality of the fundamental assumption of bimodality, sample histogram fits for three representative magnitude ranges are shown in Figure 8. Two modes are definitely required, but the arbitrary addition of more than two subcomponents results in no formal improvement to the total color distribution (see also Larsen et al. 2001; Harris 2009; Cockcroft et al. 2009, for similar tests in several other galaxies), showing that bimodality is both a necessary and sufficient condition for describing the broadband–color distribution.

Zoom In Zoom Out Reset image size

7.1. The Blue Sequence

The behavior of the key bimodal parameters (the means μ₁, μ₂ and the dispersions σ₁, σ₂) with luminosity is shown in Figure 9. A nonzero MMR will appear as a systematic trend of the mean color μ with magnitude. By using the bimodal fits to define the central colors of each mode, we also explicitly account for any overlap of the color distributions and possible asymmetric biasing due to contamination of one mode by the other. We concentrate first on the trends along the blue sequence, which to date is the only one of the two claimed to show a clear MMR from the studies of many other galaxies.

Zoom In Zoom Out Reset image size

The blue sequence in M87 shows up as a narrow, rather sharply defined feature in the CMD that is easily traced over the entire spatial region studied here. The narrow spread and large population on the blue sequence mean that the numerical solution for its position and shape is very stable. In the subsample used here for the fits, it takes up 55%-60% of the total GC population, and the p₁ ratio stays nearly constant at all magnitude levels (upper panel of Figure 10).² Comparison of the results in Tables 2 and 3 (and displayed in Figure 9) shows that no significant changes in μ₁ or σ₁ occur when the red sequence is arbitrarily required to fall at a constant color. The mean color (normalized to R = 1.0R_eff) is 〈g − i〉(blue) = 0.80 ± 0.01 over the fainter half of the i-magnitude range studied here. Its intrinsic width has a dispersion σ₁ = 0.06 ± 0.01 that remains highly consistent with magnitude, except in the faintest range i > 22.5 where it is partially broadened by photometric measurement scatter. The Maraston (2005) population models used above then show that the dispersion in color corresponds to an intrinsic rms metallicity spread σ(Fe/H) = 0.30 ± 0.08. For comparison, Harris (2009) found σ(Fe/H)(blue) = 0.27 for the average of six other BCGs, measured from the (B − I) color index.

Zoom In Zoom Out Reset image size

The bimodal fits show that the blue sequence becomes smoothly and distinctly redder with increasing luminosity, especially for i ≲ 21.5 (M_I ≲ −10). Over this range, the mean color behaves as (g − i)_norm = (0.845 ± 0.090) − (0.021 ± 0.004)(i − 20), a slope that is different from zero at the 5σ level. The resultant scaling of GC heavy-element abundance with luminosity is then Z ∼ L^{0.25 ± 0.08}. For comparison, Cockcroft et al. (2009) and Harris (2009) find an average blue-sequence scaling of Z ∼ L^{0.3 ± 0.05} for more than a dozen large E galaxies measured from a variety of both ground-based and HST data and applying to a very similar GC luminosity range, though these studies also indicate that real differences from one galaxy to the next may exist.

At fainter levels (M_I ≳ −10, corresponding roughly to 10⁶M_☉), it is not clear that any significant trend exists and that the blue sequence could better be considered vertical. In summary, these results for M87 are closely similar to those of Harris (2009) for six other BCG giants: the MMR for the blue sequence in M87 exists, but it is nonlinear in form and most strongly affects the most massive clusters. Comparison of the Megacam result with the recent thorough analysis of the deep (V − I) data from HST/ACS (Peng et al. 2009) also shows close similarities: they found a clearly nonzero trend toward redder colors over the magnitude range M_I ≲ −9 and a nearly vertical sequence below that.

The data discussed here emphasize that detecting the MMR correlation on observational grounds depends not just on using a metallicity-sensitive photometric index, but also on having a very large sample that extends to high luminosity. If we had restricted our analysis to just the inner zone (see, for example, the first panel of Figure 6), the MMR with its nonlinear form would be essentially undetectable. The Waters et al. (2009) study, which uses (V − I) colors for the clusters that fall in a single HST/ACS field, covers only the region within 2' of the center of M87. Their conclusion that no MMR is present in fact does not disagree with the present Megacam results for two reasons: first, at lower GC luminosity M_I ≳ −10 where the vast majority of the clusters fall, the MMR actually does become more nearly vertical and no trend is expected. Second, the smaller sample size in the HST/ACS data means that only a few dozen of its clusters lie above the transition level M_I ∼ −10 where the MMR starts to become significant, making it hard to distinguish above the intrinsic scatter in the sequence. A more detailed comparative analysis of the HST/ACS data is provided by Peng et al. (2009).

7.2. The Red Sequence

In Figure 11, a final version of the CMD is shown with the best-fit lines for both sequences shown. The results of this study fall into line with the most complete recent collections of data for other large galaxies (Mieske et al. 2006; Harris 2009; Cockcroft et al. 2009). Above a transition point of about 10⁶M_☉, the low-metallicity clusters become progressively more enriched at higher luminosity, scaling as Z ∼ L^0.25. This nonlinear MMR can be reasonably well understood through a basic self-enrichment model (Strader & Smith 2008; Bailin & Harris 2009).

Zoom In Zoom Out Reset image size

At the same time, a MMR of similar amplitude and reverse direction—but of much lower significance—is found along the red sequence, Z ∼ L^{−(0.3 ± 0.15)}. Interestingly, Cockcroft et al. (2009) find a mean scaling Z(red sequence) ∼L^{−(0.1 ± 0.1)} for the average of 15 giant E galaxies, a result that is only marginally different from zero but at least points in the same direction as the M87 data. Thus it may be worth considering how the red sequence could have embedded in it some form of substructure appearing as a MMR of its own in a way that does not require self-enrichment.

To discuss the possibilities, we first ask what the essential differences are between the two modes. The blue, metal-poor clusters make up a narrow, well-defined sequence that is replicated in very much the same way over and over in galaxies of all types, differing primarily just in their total numbers within different galaxies. Their mean metallicity 〈Fe/H〉 ≃ −1.4 is only weakly correlated with galaxy mass and shows very much the same internal dispersion σ(Fe/H) ≃0.3 from one system to another (e.g., Burgarella et al. 2001; Strader et al. 2004; Brodie & Strader 2006; Peng et al. 2006). In addition, the blue GCs display a consistently old ∼13 Gyr mean age with little internal age spread (e.g., Puzia et al. 2005; Brodie & Strader 2006; Beasley et al. 2008; Marin-Franch et al. 2009). Many lines of argument advanced in recent years have built up the view that the metal-poor GCs are likely to have formed in the first brief round of star formation within the "pregalactic dwarfs" that began their star formation just at the beginning of the major stage of hierarchical merging, and additional metal-poor GCs could have accreted later as part of their host dwarf satellites (see especially Searle & Zinn 1978; Harris & Pudritz 1994; Forbes et al. 1997; Côté et al. 1998; Burgarella et al. 2001; Beasley et al. 2002; Kravtsov & Gnedin 2005; Brodie & Strader 2006; Bekki et al. 2008, among many others). These small host dwarfs with their shallow potential wells and initially metal-poor gas could support only one major round of star formation and thus left behind metal-poor stars plus much unconverted gas.

By contrast, the red, metal-richer sequence is a more heterogeneous population that in some ways is harder to model in detail. Their mean metallicity 〈Fe/H〉 ≃ −0.4 correlates weakly with host galaxy luminosity (Forbes et al. 1997; Strader et al. 2004; Brodie & Strader 2006). Both their metallicity distribution and spatial distribution within the halos of their galaxies as a whole are similar to the field-star population that has a broad, metal-rich MDF and makes up the bulk of a typical giant elliptical (Geisler et al. 1996; Harris & Harris 2002; Rhode & Zepf 2004; Forte et al. 2007). That is, the metal-rich GCs can reasonably be thought of as forming during the major stage of hierarchical merging that built most of the galaxy. Hierarchical models (e.g., Beasley et al. 2002, 2003) show that multiple "in situ" starbursts, major mergers, and accretions of gas-rich satellites would all have contributed to this stage. If this basic picture is correct, then the final population of metal-rich GCs would be expected to have a wide range of ages and metallicities, and more so in bigger galaxies that experienced more starbursts and merging events. Furthermore, the red-GC populations should also differ in detail from one galaxy to another depending on their individual merger and starburst histories.

On the observational side, the evidence points increasingly to the same view that the metal-rich GC "population" is actually a composite of objects with various ages, with its subcomponent MDFs heavily overlapped and thus of nearly similar mean metallicities (e.g. Puzia et al. 2005; Kundu et al. 2005; Brodie et al. 2005; Dirsch et al. 2005; Larsen et al. 2005; Brodie & Strader 2006; Cenarro et al. 2007; Goudfrooij et al. 2007; Beasley et al. 2008; Maybhate et al. 2009). Consistent with this view, both the relative numbers and luminosity functions of the red sequence differ noticeably from one galaxy to another, so that the high-luminosity end of the red GCLF does not extend to the same maximum even in galaxies with similar total GC populations. For example, in NGC 1399 the blue and red sequences have very similar "upper ends" while in M87 the blue sequence extends much higher (Forte et al. 2007); see also the difference between the Antlia giants NGC 3258 and 3268 (Bassino et al. 2008). By contrast with most other systems, in the Hydra and Coma supergiants NGC 3311 and NGC 4874 the red sequence is the one that extends to higher luminosity, up to cluster masses ≃10⁷M_☉ and perhaps higher (Wehner et al. 2008; Harris et al. 2009). Observations of young, massive star clusters in present-day starburst galaxies (Larsen & Richtler 2000; Larsen 2009) also show that the cluster formation efficiency and the initial cluster mass function can vary with the total star formation rate.

These various pieces of evidence point to a way to understand, at least in principle, how we might obtain a red GC sequence with an observable mass/metallicity relation that could be either positive or negative. In the simplest such case, suppose that the metal-rich GCs comprise two main subpopulations that can be thought of as the products of the two biggest mergers or starbursts that assembled the main body of the galaxy. In addition, suppose that the earlier and slightly less metal-rich of these two starbursts generated its GCs with a broader mass function that extended to higher luminosity than did the later and metal-richer burst. What we would see today in the GCS as a whole would then be a red sequence with a broad MDF and a "negative MMR" at the top end. This population-based approach has nothing to do with any self-enrichment, is completely consistent with the evidence summarized above, and does not invoke any unusual physical effect. Simulated GCSs along these lines are now being constructed to test these ideas further and will be discussed in a later paper.

Wide-field photometry with CFHT/Megacam in gri has been used to study the properties of the M87 globular cluster system. Approximately 6200 GCs have been measured to a limit i ≃ 23.0 (equivalent roughly to M_I ∼ −8.5, a little less than half the total population). The primary results of this study are these:

• 1.  
  The metal-poor (blue) and metal-rich (red) cluster sequences are clearly delineated especially in the (g − i) color index. While the red sequence does not extend detectably outward to projected galactocentric distances beyond ∼50 kpc against the field contamination, the blue sequence can be seen to 100 kpc and beyond.
• 2.  
  Both the red and blue subpopulations have shallow but significant radial metallicity gradients, corresponding to heavy-element abundance scaling with distance of Z ∼ R^−0.12 (blue) and Z ∼ R^−0.17 (red).
• 3.  
  The blue sequence shows a nonlinear mass/metallicity relation: GC metallicity increases with luminosity. This relation is nearly vertical (that is, mean color is uncorrelated with metallicity) for M_I ≳ −10, while at brighter levels we find a trend roughly reproduced by Z ∼ L^0.25. This result agrees well with the findings in numerous other giant E galaxies, as well as the best current data on M87 itself. Finding the trend depends critically on being able to sample the GC population to its highest luminosity levels. The blue-sequence MMR can be physically understood as the result of GC self-enrichment during formation.
• 4.  
  The red sequence is almost twice as broad as the blue sequence in its internal scatter of colors. The formal numerical solutions show that it has a "negative MMR" at the top end that is at least as steep as the blue-sequence MMR. However, the same numerical solutions show that this trend is not strongly significant, and simply assuming a constant mean color for the red sequence does not noticeably decrease the quality of the bimodal fits.
• 5.  
  In giant galaxies like M87, the red-GC sequence is likely to be a composite population, unlike the more internally homogeneous blue sequence. Every major starburst or gas-rich merger that contributed to the formation of the main body of the galaxy should have contributed to its metal-rich GC population. If the individual major starbursts generated GCs following slightly different initial mass functions, then features such as the negative-MMR can be built in to the total red sequence that we see today.

This work was supported by the Natural Sciences and Engineering Research Council of Canada through research grants to WEH, and by the Killam Foundation of the Canada Council through a research fellowship. It is a pleasure to acknowledge the Canadian Astronomy Data Center (CADC) and the high quality of their archive products that facilitated this study.