Joel C Roediger

Laura Ferrarese

Patrick Côté

Lauren A Macarthur

Rúben Sánchez-Janssen

John P Blakeslee

Eric W Peng

Chengze Liu

Roberto Muñoz

Jean-Charles Cuillandre

Roberto Munoz

Stephen Gwyn

Simona Mei

Samuel Boissier

Alessandro Boselli

Michele Cantiello

Stéphane Courteau

Pierre-Alain Duc

Ariane Lançon

J Christopher Mihos

Thomas H Puzia

James E Taylor

Patrick R Durrell

Elisa Toloba

Puragra Guhathakurta

Hongxin Zhang

The Next Generation Virgo Cluster

Keywords: galaxies: clusters: individual (Virgo), galaxies: dwarf, galaxies: evolution, galaxies: nuclei, galaxies: star clusters: general, galaxies: stellar content

published or not. The documents may come

INTRODUCTION

Despite the complexities of structure formation in a ΛCDM Universe, galaxies are well-regulated systems. Strong evidence supporting this statement are the many fundamental relations to which galaxies adhere: star formation rate versus stellar mass or gas density (Daddi et al. 2007;Elbaz et al. 2007;Noeske et al. 2007;Kennicutt & Evans 2012), rotational velocity versus luminosity or baryonic mass for disks (Courteau et al. 2007;McGaugh 2012;Lelli et al. 2016), the fundamental plane for spheroids (Bernardi et al. 2003;Zaritsky et al. 2012), and the mass of a central compact object versus galaxy mass (Ferrarese et al. 2006;Wehner & Harris 2006; Joel.Roediger@nrc-cnrc.gc.ca Beifiori et al. 2012;Kormendy & Ho 2013), to name several. Moreover, many of these relations are preserved within galaxy groups and clusters, demonstrating that such regulation is maintained in all environments (e.g. Blanton & Moustakas 2009). This paper focusses on the relationship between color and luminosity for quiescent ["quenched"] galaxies: the so-called red sequence [RS].

First identified by de Vaucouleurs (1961) and Visvanathan & Sandage (1977), the RS represents one side of the broader phenomenon of galaxy color bimodality (Strateva et al. 2001;Blanton et al. 2003;Baldry et al. 2004;Balogh et al. 2004;Driver et al. 2006;Cassata et al. 2008;Taylor et al. 2015), the other half being the blue cloud, with the green valley separating them. Based on the idea of passively evolving stellar pop-ulations, color bimodality is widely interpreted as an evolutionary sequence where galaxies transform their cold gas into stars within the blue cloud and move to the RS after star formation ends (e.g. Faber et al. 2007). This evolution has been partly observed through the increase of mass density along the RS towards low redshift (Bell et al. 2004;Kriek et al. 2008;Pozzetti et al. 2010), although the underlying physics of quenching remains a matter of active research. The standard view of color bimodality is a bit simplistic though insofar as the evolution does not strictly proceed in one direction; a fraction of galaxies in the RS or green valley have their stellar populations temporarily rejuvenated by replenishment of their cold gas reservoirs (Schawinski et al. 2014).

Crucial to our understanding of the RS is knowing when and how it formed. The downsizing phenomenon uncovered by spectroscopic analyses of nearby early-type galaxies (ETGs; Nelan et al. 2005;Thomas et al. 2005;Choi et al. 2014) implies that the RS was built over an extended period of time [∼5 Gyr], beginning with the most massive systems (e.g. Tanaka et al. 2005). These results support the common interpretation that the slope of the RS is caused by a decline in the metallicity [foremost] and age of the constituent stellar populations towards lower galaxy masses (Kodama & Arimoto 1997;Ferreras et al. 1999;Terlevich et al. 1999;Poggianti et al. 2001;De Lucia et al. 2007). Efforts to directly detect the formation of the RS have observed color bimodality to z ∼ 2 (Bell et al. 2004;Willmer et al. 2006;Cassata et al. 2008). More recently, legacy surveys such as GOODS, COSMOS, NEWFIRM, and UltraVISTA have shown that massive quiescent galaxies [M * 3 × 10 10 M ] begin to appear as early as z = 4 (Fontana et al. 2009;Muzzin et al. 2013;Marchesini et al. 2014) and finish assembling by z = 1-2 (Ilbert et al. 2010;Brammer et al. 2011). Growth in the stellar mass density of quiescent galaxies since z = 1, on the other hand, has occured at mass scales of M * and lower (Faber et al. 2007), consistent with downsizing.

Owing to their richness, concentration, and uniform member distances, galaxy clusters are an advantageous environment for studying the RS. Moreover, their characteristically high densities likely promote quenching and therefore hasten the transition of galaxies to the RS. In terms of formation, the RS has been identified in [proto-]clusters up to z ∼ 2 (Muzzin et al. 2009;Wilson et al. 2009;Gobat et al. 2011;Spitler et al. 2012;Stanford et al. 2012;Strazzullo et al. 2013;Cerulo et al. 2016). Much of the interest in z > 0 clusters has focussed on the growth of the faint end of the RS. Whereas scant evidence has been found for evolution of either the slope or scatter of the RS (Ellis et al. 1997;Gladders et al. 1998;Stanford et al. 1998;Blakeslee et al. 2003;Holden et al. 2004;Lidman et al. 2008;Mei et al. 2009;Papovich et al. 2010, but see Hao et al. 2009 andHilton et al. 2009), several groups have claimed an elevated ratio of bright-to-faint RS galaxies in clusters up to z = 0.8, relative to local measurements (Smail et al. 1998;De Lucia et al. 2007;Stott et al. 2007;Gilbank et al. 2008;Hilton et al. 2009;Rudnick et al. 2009, see also Boselli & Gavazzi 2014 and references therein). The increase in this ratio with redshift indicates that low-mass galaxies populate the RS at later times than high-mass systems, meaning that the former, on average, take longer to quench and/or are depleted via mergers/stripping at early epochs. These results are not without controversy, however, with some arguing that the inferred evolution may be the result of selection bias, small samples, or not enough redshift baseline (Crawford et al. 2009;Lu et al. 2009;De Propris et al. 2013;Andreon et al. 2014;Romeo et al. 2015;Cerulo et al. 2016).

As a tracer of star formation activity and stellar populations, colors also are a key metric for testing galaxy formation models. Until recently, only semi-analytic models [SAMs] had sufficient statistitcs to enable meaningful comparisons to data from large surveys. Initial efforts indicated that the fraction of red galaxies was too high in models, and thus quenching too efficient, which led to suggestions that re-accretion of SNe ejecta was necessary to maintain star formation in massive galaxies (Bower et al. 2006). Since then, a persistent issue facing SAMs has been that their RSs are shallower than observed (Menci et al. 2008;González et al. 2009;Guo et al. 2011). The common explanation for this is that the stellar metallicity-luminosity relation in the models is likewise too shallow. Font et al. (2008) demonstrated that an added cause of the excessively red colors of dwarf satellites is their being too easily quenched by strangulation, referring to the stripping of halo gas. While Font et al. (2008) increased the binding energy of this gas as a remedy, Gonzalez-Perez et al. (2014) have shown that further improvements are still needed. Studies of other models have revealed similar mismatches with observations (Romeo et al. 2008;Weinmann et al. 2011), indicating that the problem is widespread.

In this paper, we use multi-band photometry from the Next Generation Virgo Cluster Survey (NGVS; Ferrarese et al. 2012) to study galaxy colors in the core of a z = 0 cluster, an environment naturally weighted to the RS. The main novelty of this work is that NGVS photometry probes mass scales from brightest cluster galaxies to Milky Way satellites (Ferrarese et al. 2016b, hereafter F16), allowing us to characterize the RS over an unprecedented factor of >10 5 in luminosity [∼10 6 M in stellar mass] and thus reach a largely unexplored part of the color-magnitude distribution [CMD]. Given the unique nature of our sample, we also take the opportunity to compare our data to galaxy formation models, which have received scant attention in the context of cluster cores.

Our work complements other NGVS studies of the galaxy population within Virgo's core. Zhu et al. (2014) jointly modelled the dynamics of field stars and globular clusters [GCs] to measure the total mass distribution of M87 to a projected radius of 180 kpc. Grossauer et al. (2015) combined dark matter simulations and the stellar mass function to extend the stellar-to-halo mass relation down to M h ∼ 10 10 M . Sánchez-Janssen et al. (2016) statistically inferred the intrinsic shapes of the faint dwarf population and compared the results to those for Local Group dwarfs and simulations of tidal stripping. Ferrarese et al. (2016a) present the deepest luminosity function to date for a rich, volume-limited sample of nearby galaxies. Lastly, Côté et al. (in preparation) and Sánchez-Janssen et al. (in preparation) study the galaxy and nuclear scaling relations, respectively, for the same sample.

In Section 2 we briefly discuss our dataset and preparations thereof. Our analysis of the RS is presented in Section 3, while Sections 4-6 focus on comparisons to previous work, compact stellar systems [CSS] and galaxy formation models. A discussion of our findings and conclusions are provided in Sections 7-8.

DATA

Our study of the RS in the core of Virgo is enabled by the NGVS (Ferrarese et al. 2012). Briefly, the NGVS is an optical imaging survey of the Virgo cluster performed with CFHT/MegaCam. Imaging was obtained in the u * g ız bands1 over a 104 deg2 footprint centered on sub-clusters A and B, reaching out to their respective virial radii (1.55 and 0.96 Mpc, respectively, for an assumed distance of 16.5 Mpc; Mei et al. 2007;Blakeslee et al. 2009). The NGVS also obtained --band imaging for an area of 3.71 deg 2 [0.3 Mpc 2 ], roughly centered on M87, the galaxy at the dynamical center of subcluster A; we refer to this as the core region. NGVS images have a uniform limiting surface brightness of ∼29 gmag arcsec -2 . Further details on the acquisition and reduction strategies for the NGVS are provided in Ferrarese et al. (2012).

This paper focuses on the core of the cluster, whose boundaries are defined as, 12 h 26 m 20 s ≤ RA (J2000) ≤ 12 h 34 m 10 s 11 • 30 22 ≤ Dec (J2000) ≤ 13 • 26 45 and encompass four MegaCam pointings [see Figure 13 of F16]. A catalog of 404 galaxies for this area, of which 154 are new detections, is published in F16, spanning the range 8.9 ≤ g ≤ 23.7 and ≥50% complete to g ∼ 22. As demonstrated there, the galaxy photometry has been thoroughly analysed and cluster membership extensively vetted for this region; below we provide a basic summary of these endeavors. A study of the CMD covering the entire survey volume will be presented in a future contribution.

Faint [g > 16] extended sources in the core were identified using a dedicated pipeline based on ring-filtering of the MegaCam stacks. Ring-filtering replaces pixels contaminated by bright, point-like sources with the median of pixels located just beyond the seeing disk. This algorithm helps overcome situations of low surface brightness sources being segmented into several parts due to contamination. The list of candidates is then culled and assigned membership probabilities by analysing SExtractor and GALFIT (Peng et al. 2002) parameters in the context of a size versus surface brightness diagram, colors and structural scaling relations, and photometric redshifts. A final visual inspection of the candidates and the stacks themselves is made to address issues of false-positives, pipeline failures, and missed detections. After this, the remaining candidates are assigned a membership flag indicating their status as either certain, likely, or possible members.

As part of their photometric analysis, F16 measured surface brightness profiles and non-parametric structural quantities in the u * g iz bands for the core galaxies with the IRAF task ELLIPSE. These data products are complemented with similar metrics from Sérsic fits to both the light profiles and image cutouts for each source [the latter achieved with GAL-FIT]. Our work is based on the growth curves deduced by applying their [non-parametric] g -band isophotal solutions to all other bands while using a common master mask. This allows us to investigate changes in the RS as a function of galactocentric radius, rather than rely on a single aperture. Driver et al. (2006) adopted a similar approach for their CMD analysis, finding that bimodality was more pronounced using core versus global colors; our results support this point [see Fig. 4]. We extract from the growth curves all ten colors covered by the NGVS, integrated within elliptical apertures having semi-major axes of a × R e,g [R e,g = g -band effective radius], where a = 0.5, 1.0, 2.0, 3.0; we also examine colors corresponding to the total light of these galaxies. Errors are estimated following Chen et al. (2010), using the magnitude differences between F16's growth curve and Sérsic analyses, and scaling values for each set of apertures by the fraction of light enclosed. These estimates should probably be regarded as lower limits, since they do not capture all sources of systematic uncertainty.

Absolute magnitudes are computed assuming a uniform distance of 16.5 Mpc (Mei et al. 2007;Blakeslee et al. 2009) for all galaxies and corrected for reddening using the York Extinction Solver 2 (McCall 2004), adopting the Schlegel et al. (1998) dust maps, Fitzpatrick (1999) extinction law, and R V = 3.07. To help gauge the intrinsic scatter along the RS, we use recovered magnitudes for ∼40k artificial galaxies injected into the image stacks [F16] to establish statistical errors in our total light measurements. A more focussed discussion of uncertainties in the NGVS galaxy photometry may be found in Ferrarese et al. (2016a) andF16.

We note that, although the NGVS is well-suited for their detection, ultra-compact dwarfs [UCDs] are excluded from our galaxy sample for two reasons. First, they have largely been omitted from previous analyses of the RS. Second, the nature of these objects is unsettled. While many are likely the remnants of tidally-stripped galaxies (e.g. Bekki et al. 2003;Drinkwater et al. 2003;Pfeffer & Baumgardt 2013;Seth et al. 2014), the contribution of large GCs to this population remains unclear. Readers interested in the photometric properties of the UCD population uncovered by the NGVS are referred to Liu et al. (2015) for those found in the core region; still, we include UCDs in our comparisions of the colors of RS galaxies and CSS in Section 5.

THE RED SEQUENCE IN THE CORE OF VIRGO

Figure 1a plots the (u * -) colors, integrated within 1.0 R e,g , of all 404 galaxies in the core of Virgo as a function of their total g -band magnitudes. One of the most striking features in this plot is the depth to which we probe galaxy colors: at its 50%-completeness limit [M g ∼ -9], the NGVS luminosity function reaches levels that have only been previously achieved in the Local Group [i.e. comparable to the Carina dSph, and only slightly brighter than Draco; Ferrarese et al. 2016a]. This is significant as integrated colors for dwarf galaxies at these scales have, until now, been highly biased to the local volume [D ≤ 4 Mpc], incomplete, and noisy (e.g. Johnson et al. 2013). The NGVS CMD therefore represents the most extensive one to date based on homogeneous photometry, spanning a factor of 2 × 10 5 in luminosity.

Also interesting about Fig. 1a is the dearth of blue galaxies in the core of Virgo. This is more apparent in Figure 1b, where we plot histograms of (u * -) in four bins of luminosity. Three of the four samples are well described as unimodal populations rather than the bimodal color distributions typically found in large galaxy surveys (e.g. Baldry et al. 2004). The absence of a strong color bimodality in Virgo's core is not surprising though (Balogh et al. 2004;Boselli et al. 2014) and suggests that most of these galaxies have been cluster members long enough to be quenched by the environment3 . The minority of blue galaxies we find may be members that are currently making their first pericentric passage or are noncore members projected along the line-of-sight. Since our interest lies in the RS, we have inspected three-color images for every galaxy and exclude from further analysis 24 that are clearly star-forming [blue points in Fig. 1a]. Also excluded are the 56 galaxies that fall below our completeness limit [grey points], 16 whose imaging suffers from significant contamination [e.g. scattered light from bright stars; green points], and 4 that are candidate remnants of tidal stripping [red points]. While we cannot rule out a contribution by reddening to the colors of the remaining galaxies, their threecolor images do not indicate a significant frequency of dust lanes.

Figure 2 plots all ten CMDs for quiescent galaxy candidates in Virgo's core, where the colors again correspond to 1.0 R e,g . Having culled the star-forming galaxies, we can straightforwardly study the shape of the RS as a function of wavelength. In each panel of Fig. 2 we observe a clear trend, whereby for M g -14, colors become bluer towards fainter magnitudes. To help trace this, we have run the Locally Weighted Scatterplot Smoothing algorithm (LOWESS; Cleveland 1979) on each CMD; these fits are represented by the red lines in the figure. The observed trends are notable given that optical colors are marginally sensitive to the metallicities of composite stellar populations with Z 0.1Z . Simple comparisons of our LOWESS curves to stellar population models suggests that, for M g -14, metallicity increases with luminosity

Galaxy mass and possible delay times likely factor into this disagreement.

along the RS [see Fig. 9]; age trends are harder to discern with the colors available to us. A metallicity-luminosity relation for RS galaxies agrees with previous work on the stellar populations of ETGs throughout Virgo (Roediger et al. 2011b) and the quiescent galaxy population at large (e.g. Choi et al. 2014). Our suggestion though is based on fairly restrictive assumptions about the star formation histories of these galaxies [i.e. exponentially-declining, starting ≥8 Gyr ago]; more robust results on age and metallicity variations along the RS in Virgo's core from a joint UV-optical-NIR analysis will be the subject of future work.

A flattening at the bright end of the RS for Virgo ETGs was first identified by Ferrarese et al. (2006) and later confirmed in several colors by Janz & Lisker (2009, hereafter JL09). This seems to be a ubiquitous feature of the quiescent galaxy population, based on CMD analyses for nearby galaxies (Baldry et al. 2004;Driver et al. 2006). This flattening may also be present in our data, beginning at M g ∼ -19, but the small number of bright galaxies in the core makes it difficult to tell. Also, this feature does not appear in colors involving the zband, but this could be explained by a plausible error in this measurement for M87 [e.g. 0.1 mag], the brightest galaxy in our sample.

The flattening seen at bright magnitudes implies that the RS is non-linear. A key insight revealed by the LOWESS fits in Fig. 2 is that the linearity of the RS also breaks down at faint magnitudes, in all colors. The sense of this non-linearity is that, for M g -14, the local slope is shallower than at brighter magnitudes, even flat in some cases [e.g. u * -g ; see Appendix]. For several colors [e.g. --ı], the LOWESS fits sug-

0.8 1.0 1.2 1.4 1.6 u * - g ′ 0.5 0.6 0.7 0.8 0.9 1.0 1.1 g ′ - i ′ 1.2 1.4 1.6 1.8 2.0 2.2 2.4 u * - r ′ 0.4 0.6 0.8 1.0 1.2 1.4 g ′ - z ′ 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 u * - i ′ 0.1 0.2 0.3 0.4 0.5 r ′ - i ′ 1.5 2.
.2 0.3 0.4 0.5 i ′ - z ′ M g ′ (mag) Figure 2.
CMDs for quiescent galaxies in Virgo's core, in all ten colors measured by the NGVS. Fluxes have been measured consistently in all five bands within apertures corresponding to 1.0 R e,g isophote of each galaxy. Black points represent individual galaxies while red lines show non-parametric fits to the data. The RS defines a color-magnitude relation in all colors that flattens at faint magnitudes, which could be explained by a constant mean age and metallicity for the lowest-mass galaxies in this region [albeit with significant scatter; but see Fig. 3]. Representative errors for the same magnitude bins as in Fig. 1 -15 since our tests did not probe brighter magnitudes. The scatter and errors, averaged within three bins of luminosity, match quite well, especially at the faintest luminosities, suggesting minimal intrinsic scatter in the colors and stellar populations of these galaxies.

gest that the behavior at the faint end of the RS may be even more complex, but the scale of such variations is well below the photometric errors [see Fig. 3]. JL09 found that the colormagnitude relation [CMR] of ETGs also changes slope in a similar manner, but at a brighter magnitude than us [M g ∼ -16.5]; we address this discrepancy in Section 4.

An implication of the faint-end flattening of the RS is that the low-mass dwarfs in Virgo's core tend to be more alike in color than galaxies of higher mass. This raises the question of whether the scatter at the faint-end of the RS reflects intrinsic color variations or just observational errors. We address this issue in Figure 3 by comparing the observed scatter in the total colors to error estimates based on the artificial galaxies mentioned in Section 2. Shown there are LOWESS fits to the data and the rms scatter about them [solid and dashed lines, respectively], and the scatter expected from photometric errors [dotted lines]. Both types of scatter have been averaged within three bins of magnitude: -15 < M g ≤ -13, -13 < M g ≤ -11, and -11 < M g ≤ -9; the comparison does not probe higher luminosities because our artificial galaxy catalog was limited to g > 16, by design. We generally find that the scatter and errors both increase towards faint magnitudes and that the two quantities match well, except in the brightest bin, where the scatter mildly exceeds the errors. For the other bins however, the intrinsic scatter must be small, strengthening the assertion that the faintest galaxies possess uniform colors [to within 0.05 mag] and, possibly, stellar populations. Deeper imaging will be needed to improve the constraints on genuine color variations at these luminosities.

The last topic examined in this section is the effect of aperture size on galaxy color. Our most important result, the flattening of the RS at faint magnitudes, is based on galaxy colors integrated within their half-light radii. Aperture effects could be significant in the presence of radial color gradients, as suggested by Driver et al. (2006), and therefore bias our inferences on the shape of the RS. In Figure 4 we show LOWESS fits to the u * -g and g -z RSs for colors measured within 0.5 R e,g , 1.0 R e,g , 2.0 R e,g , and 3.0 R e,g . These particular colors are chosen because, in the absence of deep UV and NIR photometry 4 , they provide the only leverage on stellar populations for the full NGVS dataset. We also include measurements of the scatter about these fits for the 0.5 R e,g and 3.0 R e,g apertures, represented by the shaded envelopes.

The top panel of Fig. 4 shows that u * -g changes by at most 0.04-0.06 mag at M g ≤ -17 between consecutive aperture pairs. Two-sample t-tests of linear fits to the data indicate that these differences are significant at the P = 0.01 level. Conversely, hardly any variation is seen between apertures for galaxies with M g > -16. The bottom panel of Fig. 4 demonstrates that g -z changes little with radius in most of our galaxies. Slight exceptions are the 0.5 R e,g colors for galaxies with M g ≤ -16, which differ from the 2.0 and 3.0 R e,g colors by 0.04 mag. The 1.0 R e,g colors bridge this gap, following the 0.5 R e,g sequence at M g -17 and moving towards the other sequences for brighter magnitudes.

The changes in the RS with galactocentric radius imply the existence of negative color gradients within specific regions of select galaxies. The strongest gradients are found for u *g within bright galaxies, inside 2.0 R e,g , while galaxies with M g > -15 have little-to-none in either color. Mild negative 4 UV and deep NIR imaging of the Virgo cluster exist (Boselli et al. 2011;Muñoz et al. 2014) but can only aid us for brighter galaxies and select fields, respectively. . RS in (u * -g ) and (g -z ), for different sizes of aperture used to measure galaxy colors. All four curves consider the same sample of galaxies. The choice of aperture has an impact on the slope of the RS at M g -16 mag for (u * -g ), with smaller apertures yielding steeper slopes, while the RS is more stable in (g -z ). The shaded envelopes represent the scatter about the RS for the 0.5 R e,g and 3.0 R e,g apertures. gradients are seen in g -z between 0.5 and 1.0 R e,g for galaxies with M g < -17, consistent with previous work on the spatially-resolved colors of galaxies throughout Virgo (Roediger et al. 2011a). The most important insight though from Fig. 4 is that the flattening of the RS at faint magnitudes does not apply to a specific aperture. The implications of those gradients we do detect in our galaxies, in terms of stellar populations and comparisons with galaxy formation models, will be addressed in Section 7.

COMPARISON TO PREVIOUS WORK

Before discussing the implications of our results, over the next two sections we compare our RS to earlier/ongoing work on the colors of Virgo galaxies and CSS, starting with the former. Of the several studies of the galaxy CMD in Virgo Representative errors for the NGVS are included along the bottom. (bottom) As above but restricted to the galaxies common to both samples; measurement pairs are joined with lines. The NGVS extends the CMD for this cluster faintward by ∼5 mag, with much improved photometric errors. We also find that JL09's CMR is steeper than our own at intermediate magnitudes, likely due to their inclusion of systems having recent star formation and possible errors in their sky subtraction. (Bower et al. 1992;Ferrarese et al. 2006;Chen et al. 2010;Roediger et al. 2011a;Kim et al. 2010), that of JL09 is the most appropriate for our purposes. JL09 measured colors for 468 ETGs from the Virgo Cluster Catalog (Binggeli et al. 1985), based on ugriz imaging from SDSS DR5 (Adelman- McCarthy et al. 2007). Their sample is spread throughout the cluster and has B < 18.0. Most interestingly, they showed that these galaxies trace a non-linear relation in all optical CMDs, not unlike what we find for faint members inhabiting the centralmost regions.

In Figure 5 we overlay the u * -g CMD from JL09 against our own, measured within 1.0 R e,g ; the comparison is ap-propriate since JL09 measured colors within their galaxies' r-band half-light radii. We have transformed JL09's photometry to the CFHT/MegaCam system following Equation 4 in Ferrarese et al. (2012). The top panel shows all objects from both samples, along with respective LOWESS fits, while the bottom is restricted to the 62 galaxies in common to both. We focus on the u * -g color because of its importance to stellar population analyses; indeed, this is a reason why accurate u *band photometry was a high priority for the NGVS.

The most notable feature in the top panel of Fig. 5 is the superior depth of the NGVS relative to the SDSS, an extension of ∼5 mag. There is a clear difference in scatter between the two samples, with that for JL09 increasing rapidly for M g > -18, whereas the increase for the NGVS occurs much more gradually5 [cf. Fig. 3; see Fig. 1 -16.5]. The shallower slopes found by JL09 at both ends of their CMR are seen for other colors and so cannot be explained by limitations/biases in the SDSS u-band imaging. The shallower slope at bright magnitudes substantiates what was hinted at in Fig. 2 and is more obvious in JL09 since their sample covers the full cluster6 ; the existence of this feature is also well-known from SDSS studies of the wider galaxy population (e.g. Baldry et al. 2004). The lower zeropoint of the JL09 CMR is seen in other colors too, hinting that calibration differences between SDSS DR5 and DR7 are responsible, where the NGVS is anchored to the latter (Ferrarese et al. 2012).

Lastly, the LOWESS fits in Fig. 5 indicate that, between -19 M g -16.5, the JL09 CMR has a steeper slope than the NGVS RS. This difference is significant [P = 0.01] and holds for other u * -band colors as well. This steeper slope forms part of JL09's claim that the ETG CMR flattens at M g -16.5, a feature not seen in our data. Since JL09 selected their sample based on morphology, recent star formation in dwarf galaxies could help create their steeper slope. For one, the colors of many galaxies in the JL09 sample overlap with those flagged in our sample as star-forming. Also, Kim et al. (2010) find that dS0s in Virgo follow a steeper UV CMR than dEs and have bluer UV-optical colors at a given magnitude. We therefore are unsurprised to have not observed the flattening detected by JL09.

Recent star formation cannot solely explain why JL09 find a steeper slope at intermediate magnitudes though. The bottom panel of Fig. 5 shows that, for the same galaxies, JL09 measure systematically bluer u * -g colors; moreover, this difference grows to fainter magnitudes, creating a steeper CMR. Comparisons of other colors [e.g. g -] and the agreement found therein proves that this issue only concerns JL09's u-band magnitudes. The stated trend in the color discrepancy appears inconsistent with possible errors in our SDSS-MegaCam transformations. Aperture effects can also be ruled out since the differences in size scatter about zero and never exceed 25% for any one object; besides, Fig. 4 demonstrates that color gradients in u * -g are minimal at faint magnitudes. A possible culprit may be under-subtracted backgrounds in JL09's u-band images since they performed their own sky subtraction. Therefore, we suggest that the differences between the JL09 CMR and NGVS RS for M g > -19 can be explained by: (i) a drop in the red fraction amongst Virgo ETGs between -19 M g -16.5, and (ii) JL09's measurement of systematically brighter u-band magnitudes. Despite this disagreement, these comparisons highlight two exciting aspects about the NGVS RS [and the photometry overall]: (i) it extends several magnitudes deeper than the SDSS, and (ii) the photometric errors are well-controlled up to the survey limits.

COMPARISON TO COMPACT STELLAR SYSTEMS

The NGVS is unique in that it provides photometry for complete samples of stellar systems within a single global environment, including galaxies, GCs, galactic nuclei, and UCDs. These systems are often compared to one another through their relative luminosities and sizes (e.g. Burstein et al. 1997;Misgeld & Hilker 2011;Brodie et al. 2011), whereas their relative stellar contents, based on homogeneous datasets, are poorly known. Given the depth of the NGVS RS, we have a unique opportunity to begin filling this gap by examining the colors of faint dwarfs and CSS at fixed lumuniosity.

Our samples of GCs, nuclei, and UCDs are drawn from the catalogs of Peng et al. (in preparation), F16, and Zhang et al. (2015), respectively; complete details on the selection functions for these samples may be found in those papers. Briefly though, GCs and UCDs were both identified via magnitude cuts and the u * ıK diagram (Muñoz et al. 2014), and separated from each other through size constraints [r h ≥ 11 pc for UCDs]. The validity of candidate GCs are assessed probabilistically and we use only those having a probability > 50%. All UCDs in the Zhang et al. (2015) catalog are spectroscopically-confirmed cluster members. Lastly, galactic nuclei were identified by visual inspection of the image cutouts for each galaxy and modelled in the 1D surface brightness profiles with Sérsic functions. For our purposes, we only consider those objects classified as unambiguous or possible nuclei in the F16 catalog.

In Figure 6 we plot the CMDs of galaxies and CSS in Virgo's core [left-hand side] and the color distributions for objects with g > 18 [right-hand side]; u * -g colors are shown in the upper row and g -ı in the lower. Note that we have truncated the CSS samples to 18 < g < 22 so that our comparisons focus on a common luminosity range.

An obvious difference between the distributions for galaxies and CSS at faint luminosities is the latter's extension to very red colors, whereas the former is consistent with a single color [Fig. 3]. This is interesting given that CSS have a higher surface density than the faint dwarfs in Virgo's core, suggesting that, at fixed luminosity, diffuse systems are forced to be blue while concentrated systems can have a wide range of colors. The nature of red CSS is likely explained by a higher metal content, since metallicity more strongly affects the colors of quiescent systems than age [see Fig. 9]. Also, the Spearman rank test suggests that nuclei follow CMRs in both u * -g [ρ = -0.57; p = 4 × 10 -5 ] and g -ı [ρ ∼ -0.5; p = 6 × 10 -4 ], hinting at a possible mass-metallicity relation for this population. A contribution of density to the colors of CSS is not obvious though given that many [if not most] of them were produced in the vicinity of higher-mass galaxies, and so may owe their enrichment to their local environments. The as-yet uncertain nature of UCDs as either the massive tail of the GC population or the bare nuclei of stripped galaxies also raises ambiguity on what governs their stellar contents, be it due to internal or external factors [i.e. self-enrichment versus enriched natal gas clouds].

While it is possible for CSS to be quite red for their luminosities, the majority of them have bluer colors, in both u * -g and g -ı, that agree better with those of faint RS galaxies. Closer inspection of the right-half of Fig. 6 reveals some tensions between the populations though. KS tests indicate that the null hypothesis of a common parent distribution for galaxies and GCs is strongly disfavored for u * -g and gı [p < 10 -10 ], whereas conclusions vary for UCDs and nuclei depending on the color under consideration [p u * -g ∼ 0.09 and p g -ı < 10 -4 for UCDs; p u * -g ∼ 0.007 and p g -ı ∼ 0.07 for nuclei]. The tails in the distributions for the CSS play an important role in these tests, but their removal only brings about consistency for the nuclei. For instance, clipping objects with * -g ≥ 1.2 increases the associated p-values to 0.18, 0.17, and 0.04 for UCDs, nuclei, and GCs, respectively, while p changes to ∼ 10 -4 , 0.65, and < 10 -4 by removing objects with g -ı ≥ 0.85. We have also fit skewed normal distributions to each dataset, finding consistent mean values between galaxies and CSS [except GCs, which have a larger value in g -ı], while the standard deviations for galaxies is typically larger than those for CSS. The evidence for common spectral shapes between the majority of CSS and faint galaxies in the core of Virgo is therefore conflicting. An initial assessment of the relative stellar contents within these systems, and potential trends with surface density and/or local environment, via a joint UV-optical-NIR analysis is desirable to pursue this subject further (e.g. Spengler et al., in preparation).

COMPARISON TO GALAXY FORMATION MODELS

As stated earlier, colors allow us to test our understanding of the star formation histories and chemical evolution of galaxies; scaling relations therein; and ultimately the physics governing these processes. Here we explore whether current galaxy formation models plausibly explain these subjects by reproducing the RS in the core of Virgo. The main novelty of this comparison lies in its focus on the oldest and densest part of a z ∼ 0 cluster, where members have been exposed to extreme external forces, on average, for several Gyr (Oman et al. 2013). The nature of our sample dictates that this comparison is best suited for galaxies of intermediate-to-low masses, although we still include high-mass systems for completeness. Unless otherwise stated, when discussing the slope of the RS, we are referring to the interval -19 M g -15, where its behavior is more or less linear.

We compare our results to three recent models of galaxy formation: one SAM (Henriques et al. 2015, hereafter H15) and two hydrodynamic (Illustris and EAGLE; Vogelsberger et al. 2014;Schaye et al. 2015). H15 significantly revised the L-Galaxies SAM, centered on: (i) increased efficiency of radiomode AGN feedback; (ii) delayed reincoporation of galactic winds [scaling inversely with halo mass]; (iii) reduced density threshold for star formation; (iv) AGN heating within satellites; and (v) no ram pressure stripping of hot halo gas in low-mass groups. H15 built their model on the Millenium I and II cosmological N-body simulations (Springel et al. 2005;Boylan-Kolchin et al. 2009), enabling them to produce galaxies over a mass range of 10 7 < M * < 10 12 M . Their revi- Since our intent is to compare these stellar systems within a common magnitude range, only those CSS having 18 < g < 22 are plotted. Representative errors for each population at faint magnitudes are included at bottom-left. (bottom row) As above but for the g -ı color. At faint magnitudes, comparitively red objects are only found amongst the CSS populations; their colors are likely caused by a higher metal content than those for galaxies of the same luminosity.

sions helped temper the persistent issues of SAMs having too large a blue and red fraction at high and low galaxy masses, respectively (Guo et al. 2011;Henriques et al. 2013). Illustris consists of eight cosmological N-body hydro simulations, each spanning a volume of ∼100 3 Mpc 3 , using the moving-mesh code AREPO. This model includes prescriptions for gas cooling; stochastic star formation; stellar evolution; gas recycling; chemical enrichment; [kinetic] SNe feedback; supermassive black hole [SMBH] seeding, accretion and mergers; and AGN feedback. The simulations differ in terms of the resolution and/or particle types/interactions considered; we use the one having the highest resolution and a full physics treatment. EAGLE comprises six simulations with a similar nature to Illustris but run with a modified ver-sion of the SPH code GADGET 3 instead. The simulations differ in terms of resolution, sub-grid physics, or AGN parameterization, where the latter variations produce a better match to the z ∼ 0 stellar mass function and high-mass galaxy observables, respectively. The fiducial model [which we adopt] includes radiative cooling; star formation; stellar mass loss; feedback from star formation and AGN; and accretion onto and mergers of SMBHs. Modifications were made to the implementations of stellar feedback [formerly kinetic, now thermal], gas accretion by SMBHs [angular momentum now included], and the star formation law [metallicity dependence now included]. The galaxy populations from Illustris and EA-GLE both span a range of M * 10 8.5 M .

We selected galaxies from the z = 0.0 snapshot of H15 M g (mag)

Figure 7. Comparison of the NGVS RS to those from galaxy formation models, with gray circles marking the positions of the observed galaxies. The shaded region surrounding each model curve indicates the 1-σ scatter, measured in five bins of luminosity. Curves for Illustris do not appear in panels showing u * -band colors since their subhalo catalogs lack those magnitudes. In every color, models uniformly predict a shallower slope for the RS than is observed in cluster cores.

that inhabit massive halos [M h > 10 14 M ], have non-zero stellar masses, are quenched [sSFR < 10 -11 yr -1 ] and bulgedominated [B/T > 0.5, by mass]; the last constraint aims to weed out highly-reddened spirals. We query the catalogs for both the Millenium I and II simulations, where the latter fills in the low-mass end of the galaxy mass function, making this sample of model galaxies the best match to the luminosity/mass range of our dataset. Similar selection criteria were used to obtain our samples of Illustris and EAGLE galaxies, except that involving B/T since bulge parameters are not included with either simulation's catalogs. We also imposed a resolution cut on Illustris such that each galaxy is populated by ≥240 star particles [minimum particle mass = 1.5 × 10 4 M ]. A similar cut is implicit in our EAGLE selection as SEDs are only available for galaxies having M * 10 8.5 M . Interestingly, most of the brightest cluster galaxies in EAGLE are not quenched, such that we make a second selection to incorporate them in our sample; no such issue is found with Illustris. Broadband magnitudes in the SDSS filters were obtained from all three models and transformed to the CFHT/MegaCam system [see Section 4]. We note that these magnitudes and the associated colors correspond to the total light of these galaxies.

A final note about this comparison is that we stack clusters from each model before analysing its RS. The high densities of cluster cores make them difficult to resolve within cosmological volumes, particularly for hydro simulations, leading to small samples for individual clusters. Stacking is therefore needed to enable a meaningful analysis of the model CMD for quenched cluster-core galaxies. H15, Illustris, and EAGLE respectively yield ∼15k, 144, and 157 galaxies lying within 300 kpc of their host halos' centers, which is roughly equivalent to the projected size of Virgo's core [as we define it]. Note that the much larger size of the H15 sample is explained by the greater spatial volume it models and the fainter luminosities reached by SAMs [M g ≤ -12, compared to M g -15 for hydro models].

In Figure 7 we compare the RS from Fig. 2 [black] to those from H15 [red], Illustris [green], and EAGLE [blue], where the curves for the latter were obtained in identical fashion to those for the NGVS. The shaded regions about each model RS convey the 1σ scatter within five bins of luminosity. The Illustris RS does not appear in the panels showing u * -band colors since their catalogs lack SDSS u-band magnitudes.

The clear impression from Fig. 7 is that no model reproduces the RS in Virgo's core, with model slopes being uniformly shallower than observed. Two-sample t-tests of linear fits to the data and models show that these differences are significant at the P = 0.01 level, except for the case of the EAGLE models and g -color [P = 0.09]. Further, the H15 RS exhibits no sign of the flattening we observe at faint magnitudes; the hydro models unfortunately lack the dynamic range needed to evaluate them in this regard.

The model RSs also differ from one another to varying degrees. First, H15 favors a shallower slope than the hydro models. Second, the color of the H15 RS monotonically reddens towards bright magnitudes whereas the hydro RSs turnover sharply at M g -19. EAGLE and Illustris agree well except for the ubiquitos upturn at faint magnitudes in the latter's RS [marked with dashed lines]. These upturns are created by the resolution criterion we impose on the Illustris catalog and should be disregarded. Underlying this behavior is the fact that lines of constant M * trace out an approximate anti- correlation in color-magnitude space (Roediger & Courteau 2015), a pattern clearly seen when working with richer samples from this model [e.g. galaxies from all cluster-centric radii]. Third, the scatter in H15 is typically the smallest and approximately constant with magnitude, whereas those of the hydro models are larger and increase towards faint magnitudes, more so for Illustris. Given that we find little intrinsic scatter in the NGVS RS at M g > -15 [Fig. 3], H15 appears to outperform the hydro models in this regard, although we can only trace the latter's scatter to M g ∼ -15. Other differences between Illustris and EAGLE appear for the colors g -ı, --ı, and ı-z , in terms of turnovers, slopes and/or zeropoints, all of which are significant [P = 0.01]. It is worth noting that while Fig. 7 references colors measured within 1.0 R e,g for NGVS galaxies [to maximize their numbers], the agreement is not much improved if we use colors from larger apertures. The conflicting shapes on the RS from data and models could be viewed in one of two ways: (i) the core of Virgo is somehow special, or (ii) models fail to reproduce the evolution of cluster-core galaxies. To help demonstrate that the latter is more probable, we compare the same models against a separate dataset for nearby clusters. WINGS (Fasano et al. 2002(Fasano et al. , 2006) is an all-sky survey of a complete, X-ray selected sample of 77 galaxy clusters spread over a narrow redshift range [z = 0.04 -0.07]. Valentinuzzi et al. (2011) measured the slope of the RS for 72 WINGS clusters using BV photometry for galaxies in the range -21.5 ≤ M V ≤ -18. We have done likewise for each well-populated [N > 100] model cluster, using the Blanton & Roweis (2007) filter transformations to obtain BV photometry from SDSS gr-band magnitudes.

Figure 8 compares the distribution of RS slopes from WINGS and galaxy formation models, with the dashed line in the top panel indicating the value in Virgo's core, which fits comfortably within the former. Each model distribution is shown for the two closest snapshots to the redshift limits of the WINGS sample. In the case of H15 and Illustris, these snapshots bracket the WINGS range quite well, whereas the redshift sampling of EAGLE is notably coarser. The latter fact may be important to explaining the difference between the two distributions for this model, since z = 0.1 corresponds to a look-back time of ∼1.3 Gyr. On the other hand, H15 and Illustris suggest that the RS slope does not evolve between z = 0.07/0.08 and 0.03. We have not tried to link model clusters across redshifts as parsing merger trees lies beyond the scope of this work. Observations though support the idea of a static slope in clusters over the range z = 0 -1 (Gladders et al. 1998;Stanford et al. 1998;Blakeslee et al. 2003;Ascaso et al. 2008).

Fig. 8 demonstrates that the distributions for the WINGS and model clusters are clearly incompatible, with the models, on average, preferring a shallower slope for the RS. The sense of this discrepancy is the same as that seen in Fig. 7 between the core of Virgo and the models. A caveat with the comparisons to WINGS though is that the model slopes have all been measured in the respective rest-frames of the clusters. In other words, the model slopes could be biased by differential redshifting of galaxy colors as a function of magnitude [e.g. fainter galaxies reddened more than brighter ones]. To address this, we have simulated the effect of k-corrections using the median of the EAGLE distribution at z = 0.1, finding it would steepen this cluster's RS by -0.01. While significant, we recall that the redshift range for the WINGS sample is z = 0.04 -0.07, such that the mean k-correction to the model slopes is likely smaller than this value and would therefore not bring about better agreement.

Given the value of the above comparisons for testing galaxy formation models, we provide in the Appendix parametric fits to the NGVS RS in every color [measured at 1 R e,g ]. These fits reproduce our LOWESS curves well and enable the wider community to perform their own comparisons.

DISCUSSION

Figure 1 indicates that >90% of the galaxy population within the innermost ∼300 kpc of the Virgo cluster has likely been quenched of star formation. This makes the population ideal for studying the characteristics of the RS, such as its shape and intrinsic scatter. Our analysis demonstrates that, in all optical colors, the RS is (a) non-linear and (b) strongly flattens in the domain of faint dwarfs. The former behavior had already been uncovered in Virgo, albeit at the bright end (Ferrarese et al. 2006;JL09), while the latter, which is new, begins at -14 < M g < -13 [see Appendix], well above the completeness limit of the NGVS. No correlation is observed between color and surface brightness, in bins of luminosity, for M g > -15, implying that the faint-end flattening is not the result of bias or selection effect.

The RS follows the same general shape at M g < -14 in each color, which may have implications for trends in the stellar populations of these galaxies. Assuming that bluer [e.g. u * -g ] and redder [e.g. g -z ] colors preferentially trace mean age and metallicity (Roediger et al. 2011b), respectively, the decrease in color towards faint magnitudes over the range -19 M g ≤ -14 hints that the populations become younger and less enriched (consistent with downsizing; Nelan et al. 2005), with two exceptions. The flattening at bright magnitudes, seen better in samples that span the full cluster (JL09) and the global galaxy population (Baldry et al. 2004), signals either a recent burst of star formation within these galaxies or an upper limit to galactic chemical enrichment. The latter seems more likely given that the stellar mass-metallicity relation for galaxies plateaus at M * 10 11.5 M (Gallazzi et al. 2005). The other exception concerns the flattening at the faint-end of the RS.

7.1. What Causes the Faint-End Flattening of the RS? If colors reasonably trace stellar population parameters [see next sub-section], then arguably the most exciting interpretation suggested by the data is that the faint dwarfs in Virgo's core have a near-uniform age and metallicity, over a range of ∼3-4 magnitudes. This would imply that the known stellar population scaling relations for quiescent galaxies of intermediate-to-high mass (e.g. Choi et al. 2014) break down at low masses [below ∼4 × 10 7 M ; see Appendix] and, more fundamentally, that the physics governing the star formation histories and chemical enrichment of galaxies decouples from mass at these scales.

Given the nature of our sample, the above scenario begs the questions of whether the faint-end flattening of the RS is caused by the environment, and if so, when and where the quenching occurs. While Geha et al. (2012) make the case that dwarfs with M * < 10 9 M must essentially be satellites in order to quench (also see Slater & Bell 2014;Phillips et al. 2015;Davies et al. 2016), we know little of the efficiency and timescale of quenching at low satellite masses and as a function of host halo mass. Using Illustris, Mistani et al. (2016) showed that, on average, the time to quench in low-mass clusters decreases towards low satellite masses, from ∼5.5 Gyr to ∼3 Gyr, over the range 8.5 log M * 10. Slater & Bell (2014) combine measurements of Local Group dwarfs with N-body simulations to suggest that, in such groups, galaxies of M * 10 7 M quench within 1-2 Gyr of their first pericenter passage. However, Weisz et al. (2015) compared HST/WFPC2 star formation histories to predicted infall times based on Via Lactea II (Diemand et al. 2008), finding that many dwarfs in the Local Group likely quenched prior to infall.

In addition to reionization, pre-processing within smaller host halos may play a key role in explaining why many Local Group dwarfs ceased forming stars before their accretion. Likewise, pre-processing must also be considered when trying to understand issues pertaining to quenching of cluster galaxies (e.g. McGee et al. 2009;De Lucia et al. 2012;Wetzel et al. 2013;Hou et al. 2014;Taranu et al. 2014), such as the cause of Virgo's flattened RS at faint magnitudes. Wetzel et al. (2013) deduced where satellites of z = 0 groups/clusters were when they quenched their star formation, by modelling SDSS observations of quiescent fractions with mock catalogs. They found that for host halo masses of 10 14-15 M the fraction of satellites that quenched via pre-processing increases towards lower satellite masses, down to their completeness limit of M * ∼ 7 × 10 9 M , largely at the expense of quenching in-situ. Extrapolating this trend to lower satellite masses suggests that the majority of the quiescent, low-mass dwarfs in Virgo were quenched elsewhere. This suggestion is consistent with abundance matching results for our sample (Grossauer et al. 2015), which indicate that only half of the core galaxies with M * = 10 6-7 M were accreted by z ∼ 1 (see also Oman et al. 2013).

Assuming that the flattening of the RS reflects an approximate homogeneity in stellar contents [i.e. constant mean age] and isolated low-mass dwarfs have constant star formation histories (e.g. Weisz et al. 2014), then the low-mass dwarfs in Virgo's core must have quenched their star formation coevally. Moreover, when coupled with a significant contribu- tion by pre-processing, it is implied that these galaxies are highly susceptible to environmental forces, over a range of host masses. This seems plausible given the very high quiescent fractions [>80%] for satellites between 10 6 < M * /M < 10 8 within the Local Volume (Phillips et al. 2015), which has led to the idea of a threshold satellite mass for effective environmental quenching (Geha et al. 2012;Slater & Bell 2014).

If synchronized quenching of low-mass dwarfs in groups [at least to ∼10 12 M ] leads to a flattened faint-end slope of the core RS, we should expect to find the same feature for dwarfs drawn from larger cluster-centric radii. This follows from the fact that a satellite's cluster-centric radius correlates with its infall time (De Lucia et al. 2012) and that the fraction of satellites accreted via groups increases towards low redshift (McGee et al. 2009). Studying the properties of the RS as a function of cluster-centric position (e.g. see Sánchez-Janssen et al. 2008) will be the focus of a future paper in the NGVS series.

Caveats

A major caveat with the above interpretations is that optical colors are not unambiguous tracers of population parameters, especially at low metallicities (Conroy & Gunn 2010). To this point, Kirby et al. (2013) have shown that stellar metallicity increases monotonically for galaxies from [Fe/H] ∼ -2.3 at M * = 10 4 M to slightly super-solar at M * = 10 12 M . Assuming this trend holds in all environments, we can check for any conditions under which the RS would flatten at faint magnitudes. In the middle and bottom panels of Figure 9 we compare the u * -g and g -z color-mass relations in Virgo's core [black lines] to those predicted by the Flexible Stellar Population Synthesis [FSPS] model (Conroy et al. 2009), where the Kirby et al. relation [top panel] is used to assign masses to each model metallicity track and lines of constant age are colored from purple [∼2 Gyr] to red [∼15 Gyr]. Other models (e.g. Bruzual & Charlot 2003) prove inadequate for our needs due to their coarse sampling of metallicity space over the range Z ∼ 4 × 10 -4 to 4 × 10 -3 . Error bars on the NGVS relations reflect standard errors in the mean, measured within seven bins of luminosity [having sizes of 0.5-2.0 dex]. Although we assume single-burst star formation histories for this test, qualitatively similar trends are expected for more complex treatments (e.g. constant star formation with variable quenching epochs; Roediger et al. 2011b).

Since the intent of Fig. 9 is to explore an alternative interpretation of the faint-end flattening of the RS, we limit our discussion to the range M * < 10 8 M , but show the full relations for completeness. Within that range, we find that the data are indeed consistent with Kirby et al.'s mass-metallicity relation, provided that age does not vary greatly therein. Moreover, the color-mass relation for select ages transitions to a flatter slope at lower masses. This confirms our previous statement that it is difficult to meaningfully constrain metallicities below a certain level with optical colors [Z 10 -3 in the case of FSPS], even when ages are independently known. The inconsistent ages we would infer from the the u * -g and g -z colors could likely be ameliorated by lowering the zeropoint of the Kirby et al. relation since the former color responds more strongly to metallicity for log(Z/Z ) -1. The comparisons shown in Fig. 9 therefore cast doubt on whether the flattening of the RS at faint magnitudes implies both a constant age and metallicity for cluster galaxies at low masses. Distinguishing between these scenarios will be more rigorously addressed in forth-coming work on the stellar populations of NGVS galaxies that incorporates UV and NIR photometry as well.

Shortcomings of Galaxy Formation Models

Regardless of the uncertainties inherent to the interpretation of optical colors, we should expect galaxy formation models to reproduce our observations if their physical recipes are correct. Our test of such models is special in that it focuses on the core of a z = 0 galaxy cluster, where the time-integrated effect of environment on galaxy evolution should be maximal. However, Fig. 7 shows that current models produce a shallower RS than observed, in all colors. This issue is not limited to Virgo's core, as Fig. 8 demonstrates that the distributions of RS slopes for entire model clusters populate shallower values than those measured for other nearby clusters. On a related note, Licitra et al. (2016) have shown that clusters at z < 1 in SAMs suffer from ETG populations with too low an abundance and too blue colors, while ∼10% of model clusters have positive RS slopes. On the other hand, Merson et al. (2016) found broad consistency between observations and SAMs in the zeropoint and slope of the RS in z > 1 clusters. This suggests that errors creep into the evolution of cluster galaxies in SAMs at z < 1.

The discrepancies indicated here follow upon similar issues highlighted by modellers themselves. H15 showed that their model produces a RS having bluer colors than observed in the SDSS for galaxies with M * ≥ 10 9.5 M . Vogelsberger et al. (2014) found the Illustris RS suffers the same problem, albeit at higher masses [M * > 10 10.5 M ], while also producing too low of a red fraction at M * < 10 11 M . Trayford et al. (2015) analyzed the colors of EAGLE galaxies, finding that its RS matches that from the GAMA survey (Taylor et al. 2015) for M r < -20.5, but is too red at fainter magnitudes. Our comparisons build on this work by drawing attention to model treatments of dense environments over cosmic time and [hopefully] incentivize modellers to employ our dataset in future work, especially as they extend their focus towards lower galaxy masses. To this end, the reader is reminded of the parametric fits to the NGVS RS provided in the Appendix.

Naturally, the root of the above discrepancies is tied to errors in the stellar populations of model galaxies. The supplementary material of H15 shows that the model exceeds the mean stellar metallicity of galaxies over the range 10 9.5 < M * 10 10 M by several tenths of a dex while undershooting measurements at 10 10.5 < M * 10 11 M by ∼0.1-0.2 dex. The issues with the H15 RS then seems to reflect shortcomings in both the star formation and chemical enrichment histories of their model galaxies. Part of the disagreement facing Illustris likely stems from the fact that their galaxies have older stellar populations than observed, by as much as 4 Gyr, for M * 10 10.5 M (Vogelsberger et al. 2014). Schaye et al. (2015) showed that EAGLE produces a flatter stellar mass-metallicity relation than measured from local galaxies due to too much enrichment at M * 10 10 M . Our inspection of the stellar populations in H15 and EAGLE reveals that their cluster-core galaxies, on average, have roughly a constant mass-weighted age [∼10-11 Gyr] and follow a shallow mass-metallicity relation, with EAGLE metallicities exceeding H15 values by ∼0.3 dex 7 . The discrepant colors produced by models thus reflect errors in both the star formation histories and chemical enrichment of cluster galaxies; for instance, 7 We omit Illustris from this dicussion as their catalogs do not provide mean stellar ages of their galaxies.

ram pressure stripping may be too effective in quenching cluster dwarfs of star formation (e.g. Steinhauser et al. 2016).

Two critical aspects of the RS that modellers must aim to reproduce are the flattenings at both bright and faint magnitudes. The former is already a contentious point between models themselves, with hydro varieties producing a turnover and while SAMs continuously increase [Fig. 7]. We remind the reader that our LOWESS curves are too steep for M g -19 since they essentially represent an extrapolation from intermediate magnitudes; the bright-end flattening is clearly visible in other datasets that span the full cluster and contain more of such galaxies [Fig. 5]. Hydro models appear to supercede SAMs in this regard, although it may be argued that their turnovers are too sharp. In the case of EAGLE, however, it is unclear what causes this turnover as several of their brightest cluster galaxies are star-forming at z = 0 while their luminosity-metallicity relation inverts for M g ≤ -20.

At present, only SAMs have the requisite depth to check for the flattening seen at the faint end of the RS; the effective resolution of cosmological hydro models is too low to probe the luminosity function to M g ∼ -13. Fig. 7 shows that the H15 RS exhibits no obvious change in slope at faint magnitudes, let alone the pronounced flattening seen in Virgo. The faintend flattening is a tantalizing feature of the RS that may hold new physical insights into the evolution of cluster galaxies of low mass. Addressing the absence of these features should be a focal point for future refinements of galaxy formation models.

CONCLUSIONS

We have used homogeneous isophotal photometry in the u * g iz bands for 404 galaxies belonging to the innermost ∼300 kpc of the Virgo cluster to study the CMD in a dense environment at z = 0, down to stellar masses of ∼ 10 6 M . Our main results are:

• The majority of galaxies in Virgo's core populate the RS [red fraction ∼ 0.9];

• The RS has a non-zero slope at intermediate magnitudes [-19 < M g < -14] in all colors, suggesting that stellar age and metallicity both decrease towards lower galaxy masses, and has minimal intrinsic scatter at the faint end;

• The RS flattens at both the brightest and faintest magnitudes [M g < -19 and M g > -14, respectively], where the latter has not been seen before;

• Galaxy formation models produce a shallower RS than observed at intermediate magnitudes, for both Virgo and other nearby clusters. Also, the RS in hydrodynamic models flattens for bright galaxies while that in SAMs varies monotonically over the full range of our dataset.

The flattening of the RS at faint magnitudes raises intriguing possibilities regarding galaxy evolution and/or cluster formation. However, these hinge on whether the flattening genuinely reflects a homogeneity of stellar populations in low-mass galaxies or colors becoming a poor tracer of age/metallicity at low metallicities [e.g. log(Z/Z ) -1.3]. This issue will be addressed in a forthcoming paper on the stellar populations of NGVS galaxies. A topic worth exploring with our parametric fits is whether the flattening of the RS occurs at a common magnitude for all colors. This can be done with the parameter M g ,0 and Table 1 shows that -14 ≤ M g ,0 ≤ -13, in a large majority of cases. For g -˚and ı-z the transition magnitude is brighter than -14, which might be explained by the fact that these colors sample short wavelength baselines and that the RS spans small ranges therein [∼0.25 and 0.15 mag, respectively]. It is also likely that the posterior distributions for the parameters in our fit are correlated.

Another way to assess the magnitude at which the RS flattens involves measuring the local gradient along our LOWESS fits, the advantage being that this approach is non-parametric. Figure 11 shows our RSs [top panel], scaled to a common zeropoint [arbitrarily established at M g ∼ -14], and the variation of the local gradient as a function of magnitude [bottom panel]. We measure the local gradient using a running bin of 9 [thin line] or 51 [thick line] data points, with the smaller bin allowing us to extend our measurements to brighter magnitudes, where our sample is sparse.

The local gradient varies in a consistent way for all colors at M g ≤ -12: the gradient is roughly constant and negative at bright magnitudes and becomes more positive towards faint magnitudes. The behaviors of the gradients at M g > -12 are more irregular as small fluctuations in the LOWESS curves are amplified when the gradients hover near zero. These behaviors are beyond this discussion however; we are interested in the locations where the rate of change of the gradients is maximized [i.e. the second derivatives of the RSs peak]. Disregarding the curves at M g > -12 then, the bottom panel of Fig. 11 shows that the rate of change maximizes in the range -14 < M g < -13, corresponding to an approximate stellar mass of ∼ 4 × 10 7 M (Ferrarese et al. 2016a). The approximate synchronicity of the flattening of the RS adds further insight to our main result on the flattening of the RS by suggesting a mass scale below which internal processes may cease to govern the stellar populations and evolution of dwarf satellites.

Figure 1 .

 1 Figure 1. (a) (u * -) color versus absolute g -band magnitude for the 404 galaxies in the core of Virgo. Colored points are purged from our sample of RS candidates due to obvious star formation activity [blue], our completeness limits [grey], significant image contamination [green], or suspected tidal stripping [red]. The vertical lines indicate bins of magnitude referenced in the right-hand panel, with representative errors plotted in each. (b) Color distributions within the four magnitude bins marked at left. The NGVS photometry enables a deep study of the galaxy CMD and we verify that the core of Virgo is highly deficient in star-forming galaxies.

  Figure2. CMDs for quiescent galaxies in Virgo's core, in all ten colors measured by the NGVS. Fluxes have been measured consistently in all five bands within apertures corresponding to 1.0 R e,g isophote of each galaxy. Black points represent individual galaxies while red lines show non-parametric fits to the data. The RS defines a color-magnitude relation in all colors that flattens at faint magnitudes, which could be explained by a constant mean age and metallicity for the lowest-mass galaxies in this region [albeit with significant scatter; but see Fig.3]. Representative errors for the same magnitude bins as in Fig.1are shown in each panel.

Figure 3 .

 3 Figure 3. Comparison of the observed scatter [dashed lines] about the RS [solid lines] to photometric errors [dotted lines] established from artificial galaxy tests. The comparision is limited to M g-15 since our tests did not probe brighter magnitudes. The scatter and errors, averaged within three bins of luminosity, match quite well, especially at the faintest luminosities, suggesting minimal intrinsic scatter in the colors and stellar populations of these galaxies.

  Figure4. RS in (u * -g ) and (g -z ), for different sizes of aperture used to measure galaxy colors. All four curves consider the same sample of galaxies. The choice of aperture has an impact on the slope of the RS at M g -16 mag for (u * -g ), with smaller apertures yielding steeper slopes, while the RS is more stable in (g -z ). The shaded envelopes represent the scatter about the RS for the 0.5 R e,g and 3.0 R e,g apertures.

Figure 5 .

 5 Figure 5. (top) Comparison of the u * -g CMD from JL09 for Virgo ETGs [black circles] to that measured here [red dots]. The full sample is plotted for each dataset and LOWESS fits for both are overlaid [solid lines].Representative errors for the NGVS are included along the bottom. (bottom) As above but restricted to the galaxies common to both samples; measurement pairs are joined with lines. The NGVS extends the CMD for this cluster faintward by ∼5 mag, with much improved photometric errors. We also find that JL09's CMR is steeper than our own at intermediate magnitudes, likely due to their inclusion of systems having recent star formation and possible errors in their sky subtraction.

Figure 6 .

 6 Figure6. (top row) u * -g CMD and color distributions for galaxies [circles], GCs [dots], UCDs [squares], and galactic nuclei [diamonds] within the core of Virgo. Since our intent is to compare these stellar systems within a common magnitude range, only those CSS having 18 < g < 22 are plotted. Representative errors for each population at faint magnitudes are included at bottom-left. (bottom row) As above but for the g -ı color. At faint magnitudes, comparitively red objects are only found amongst the CSS populations; their colors are likely caused by a higher metal content than those for galaxies of the same luminosity.

Figure 8 .

 8 Figure 8. Comparison of RS slopes in real (top panel) and model clusters (other panels). The model slopes are measured from those snapshots which most closely bracket the redshift range of the WINGS clusters [0.03 ≤ z ≤ 0.07]. In all cases the typical slope within model clusters is shallower than observed. The dashed line indicates the RS in Virgo's core.

Figure 9

 9 Figure 9. u * -g and g -z color-mass relations [middle and bottom panels; black lines] versus those predicted by the FSPS stellar population model [colored lines], constrained by the Kirby et al. (2013) mass-metallicity relation [top panel]. Each model relation corresponds to a certain fixed age, ranging between ∼2 Gyr [purple] and ∼15 Gyr [red] in steps of 0.025 dex. Error bars on the NGVS relations represent standard errors in the mean within bins of luminosity.

Figure 10 .

 10 Figure10. Parameteric fits [green lines] to the RS in Virgo's core, corresponding to the 1.0 R e,g -colors of NGVS galaxies. These fits are compared to the data themselves [black points] as well as non-parametric [LOWESS] fits. Points clipped from the each fit are shown in blue.

Figure 11 .

 11 Figure 11. (top) LOWESS fits from Fig. 10, scaled to a common zeropoint at M g ∼ -14. (bottom) Local gradient measured along each RS shown in the top panel using a rolling bin of either 9 [thin line] or 51 [thick line] data points; the former bin size allows us to extend our measurements up to bright galaxies. In all cases, the local gradient begins to flatten in the vicinity of M g ∼ -15.

  of Ferrarese et al. 2016a as well]. Furthermore, the JL09 CMR has a lower zeropoint [by ∼0.06] and a shallower slope than the NGVS RS for M g -19, which two-sample t-tests verify as significant [P = 0.01]. The JL09 data also exhibit a flattening of the CMR in the dwarf regime, but at a brighter magnitude than that seen in ours [M g ∼

Table 1

 1 Parameters of double power-law fit to the NGVS RS.

	Color	M g ,0	C 0	β 1	β 2	α	rms
		(mag) (mag)				(mag)
	(1)	(2)	(3)	(4)	(5)	(6)	(7)
	u * -g u * -˚-13.62 1.552 3.871 0.000 11.51 0.091 -13.52 1.040 2.624 0.000 15.98 0.078 u * -ı -13.45 1.787 4.577 0.000 11.81 0.116
	u * -z	-13.95 1.927 5.494 0.773 20.73 0.157
	g -˚-14.57 0.522 1.036 0.392 57.14 0.047
	g -ı	-13.81 0.751 1.685 0.578 1333. 0.058
	g -z	-13.74 0.852 2.808 0.413 23.39 0.124
	--ı	-13.07 0.230 0.735 0.130 96.97 0.050
	--z	-13.40 0.342 1.851 0.000 11.86 0.108
	ı-z	-14.15 0.107 1.133 0.000 15.92 0.102
	shape of the RS.

Note that the filters used in the NGVS are not identical to those of the Sloan Digital Sky Survey (SDSS;York et al.

2000), with the u * -band being the most different. Unless otherwise stated, magnitudes are expressed in the MegaCam system throughout this paper.

http://www4.cadc-ccda.hia-iha.nrc-cnrc.gc.ca/community/YorkExtinctionSolver/

The timescale associated with environmental quenching appears contentious, with some groups favoring shorter values (<2 Gyr;Boselli & Gavazzi 2014, and references therein;Haines et al. 2015) and others longer (several Gyr; e.g.Balogh et al. 2000;De Lucia et al. 2012;Taranu et al. 2014).

While the scatter in the JL09 data is likely dominated by the shallower depth of the SDSS imaging, a contribution by distance uncertainties cannot be ruled out, since the Virgo Cluster Catalog spans several sub-groups whose relative distances can exceed 10 Mpc(Mei et al. 2007).

Virgo comprises two large sub-clusters and several massive groups, such that its bright galaxies are spread throughout the cluster.

APPENDIX

Here we present parametric fits for the RS in Virgo's core based on the colors of our galaxies within 1.0 R e,g . Our purpose is to enable the wider community, particularly modellers, to compare our results to their own through simple [continuous] fitting functions. Motivated by the non-parameteric fits in Fig. 2, we choose a double power-law to describe the shape of the RS; we acknowledge that this choice is made strictly on a phenomenological basis and lacks physical motivation. This function is parameterized as,

where β 1 and β 2 represent the asymptotic slopes towards bright and faint magnitudes, respectively, while M g ,0 and C 0 correspond to the magnitude and color of the transition point between the two power-laws, and α reflects the sharpness of the transition.

We fit Equation 1 to our data through an iterative non-linear optimization of χ 2 following the L-BFGS-B algorithm (Byrd et al. 1995;Zhu et al. 1997), restricting α, β 1 , and β 2 to positive values, and M g ,0 and C 0 to lie in the respective ranges [-20, -8] and [0,20]. At each iteration, > 3σ outliers are clipped from each CMD; doing so allows the fits to better reproduce our LOWESS curves. We generally achieve convergence after 5-6 iterations while the fraction of clipped points is <10% in all cases.

Our power-law fits [green curves] are compared to the data [black points] and LOWESS fits [red curves] in Figure 10, while clipped data are represented by the blue points. The best-fit parameters are summarized in Table 1, where the final column lists the rms of each fit. Inspection of the rms values and the curves themselves indicates that our parametric fits do well in tracing the