NGVS-IR, I. A New Near-UV/Optical/Near-IR Globular Cluster Selection Tool

# The Next Generation Virgo Cluster Survey – Infrared (NGVS-IR)1: I. A new Near-UV/Optical/Near-IR Globular Cluster Selection Tool

## Abstract

The NGVS-IR project (Next Generation Virgo Survey – Infrared) is a contiguous near-infrared imaging survey of the Virgo cluster of galaxies. It complements the optical wide-field survey of Virgo (NGVS). The current state of NGVS-IR consists of -band imaging of 4 centered on M87, and and -band imaging of covering the region between M49 and M87. In this paper, we present the observations of the central 4 centered on Virgo’s core region. The data were acquired with WIRCam on the Canada-France-Hawaii Telescope and the total integration time was 41 hours distributed in 34 contiguous tiles. A survey-specific strategy was designed to account for extended galaxies while still measuring accurate sky brightness within the survey area. The average limiting magnitude is  AB mag and the 50% completeness limit is  AB mag for point source detections, when using only images with better than 0.7″ seeing (median seeing ). Star clusters are marginally resolved in these image stacks, and Virgo galaxies with  AB mag arcsec are detected. Combining the data with optical and ultraviolet data, we build the color-color diagram which allows a very clean color-based selection of globular clusters in Virgo. This diagnostic plot will provide reliable globular cluster candidates for spectroscopic follow-up campaigns needed to continue the exploration of Virgo’s photometric and kinematic sub-structures, and will help the design of future searches for globular clusters in extragalactic systems. We show that the new diagram displays significantly clearer substructure in the distribution of stars, globular clusters, and galaxies than the diagram – the NGVS+NGVS-IR-equivalent of the diagram, which is widely used in cosmological surveys. Equipped with this powerful new tool, future NGVS-IR investigations based on the diagram will address the mapping and analysis of extended structures and compact stellar systems in and around Virgo galaxies.

galaxies: clusters: individual (Virgo) – galaxies: distances and redshifts – galaxies: luminosity function, mass function – galaxies: photometry – galaxies: star clusters: general
2

## 1. Introduction

The proximity of the Virgo galaxy cluster, located at a distance of 16.50.2 Mpc (Mei et al., 2007; Blakeslee et al., 2009), makes it one of the closest and most thoroughly studied laboratories of star and galaxy formation in the nearby universe. With about 2000 known galaxy members (Binggeli et al., 1985; Gavazzi et al., 2003), more than globular clusters (GCs) in M 87 and M 49 alone (Tamura et al., 2006a, b; Peng et al., 2008) and a total GC population of for the entire NGVS region (Durrell et al. 2013, in prep.), it offers vast sample statistics to investigate the evolution of stellar populations and the formation and assembly histories of their host galaxies in a relatively dense cluster environment.

The first detailed census of the Virgo galaxy cluster was published by Ames (1930) and consisted of a catalog of 2278 galaxies brighter than 18th magnitude. This survey led to the identification of several background clusters and to the discovery of a Southern extension of the Virgo cluster, which are impressive results given the technical limitations existing at that time. A quarter of a century later, Reaves (1956) published the first comprehensive catalog of dwarf galaxies in the Virgo cluster based on the study of photographic plates taken with the Lick 20-inch astrograph. This work showed the presence of a large number of dwarf elliptical galaxies beyond the local group and started a new phase of the Virgo cluster research. The most important contribution to assessing Virgo cluster membership was published by Binggeli et al. (1985) in the form of the famous Virgo Cluster Catalog (VCC). This catalog consists of 2096 galaxies, covers an area of and is complete down to mag.

### 1.1. The Next Generation Virgo Cluster Survey

The VCC has been the reference standard of Virgo for more than a quarter of a century and it is about to be superseded by the Next Generation Virgo Cluster Survey (NGVS, Ferrarese et al., 2012). The NGVS is a multi-passband optical survey conducted with MegaCam at the Canada-France-Hawaii Telescope (CFHT) on Mauna Kea, Hawaii, a 3.6-meter telescope that benefits from superb seeing and transparency conditions.

The NGVS survey area covers 104 contiguous square degrees out to the virial radius around M 87 and M 49 to unprecedented photometric depths in the five optical filters3 and . With substantially sub-arcsecond seeing, a surface brightness sensitivity better than  mag arcsec and a point-source detection limit of  mag (), this survey will extend our optical view of the galaxy luminosity function and the scaling relations of galaxies to luminosities more than five magnitudes below the faintest VCC galaxies4.

The five optical passbands of NGVS combined with the superior spatial resolution of the images (-band seeing is better than  pc at Virgo distance) provide a first classification of the sources, but leave important areas of ambiguity. In particular, the identification and study of compact stellar systems and their stellar population properties is challenged by the age-metallicity-extinction degeneracy of optical colors (Worthey, 1994; Arimoto, 1996). The difficulties begin already at the level of sample definition: The distribution of old metal-poor GCs merges into the parameter space of dwarf stars in optical color-color diagrams, while metal-rich GCs overlap with the locus of remote late-type galaxies. Around selected Virgo galaxies, the pointed observations of the HST/ACS Virgo Cluster Survey (ACSVCS) have allowed a less ambiguous determination of the GC population due to the improved spatial resolution of HST (Côté et al., 2004; Jordán et al., 2009). Spectroscopic surveys are progressively extending the collection of well-classified objects to vaster areas, though with inevitable sparseness and depth limitations (e.g. Park et al., 2012; Romanowsky et al., 2012; Schuberth et al., 2012; Pota et al., 2013).

### 1.2. The Near-Infrared Complement: NGVS-IR

The NGVS near-infrared project (NGVS-IR) is designed to complement the NGVS with deep, spatially well-resolved images at wavelengths longer than the band to provide wider sampling of the spectral energy distribution (SED). The first part of this project, and the subject of this paper, is a deep -band survey conducted with the Wide-field InfraRed Camera (WIRCam; Puget et al., 2004) at CFHT that covers the central square degrees around the Virgo cluster core, the location of the cD galaxy M 87 and its vicinity. This area was the first fully surveyed in as part of the NGVS pilot survey, and we will refer to this field as the pilot field hereafter. The NGVS-IR observations are currently being expanded beyond the pilot field with the 4-meter class Visible and Infrared Survey Telescope for Astronomy (VISTA; Dalton et al., 2006) at the European Southern Observatory (ESO) on Cerro Paranal in Chile. These observations slightly overlap with the pilot field and cover regions around the brightest elliptical galaxy in Virgo, M 49 in the filters and , as well as a strip of and -band images that connects the pilot field and the M 49 field, and augments the coverage of the M 86 group. Taken altogether, the current NGVS-IR survey area provides a total contiguous spatial coverage of the central regions in Virgo of more than deg.

### 1.3. The Need for NGVS-IR

The intermediate goals of the NGVS-IR survey are the development of a photometric selection method to identify globular clusters (GCs) in the Virgo core region, the study of the age-metallicity-mass relations of GCs in M87 and the intracluster medium, the characterization of the stellar populations of Ultra-Compact Dwarf galaxies (UCDs) by combining optical/near-IR colors, the measurement of accurate distances to bright Virgo galaxies by using the Surface Brightness Fluctuation (SBF) method (Tonry et al., 1990), the study of the star formation histories and mass profiles of several Virgo galaxies, and lastly, the determination of the galaxy mass function of the Virgo cluster.

Although the Virgo cluster has already been observed at near-IR wavelengths, none of the available near-IR surveys is deep enough to accomplish the NGVS science goals. The Two Micro All Sky Survey (2MASS, Skrutskie et al., 2006) provides uniform, precise photometry and astrometry over the entire celestial sphere in and bandpasses with limiting magnitudes of 15.8, 15.1 and 14.3 mag, respectively. The 2MASS survey offers the largest coverage of the Virgo cluster in near-IR wavelengths, but its shallow photometry is not well suited for deep extragalactic studies. Deeper photometry is provided by the UKIRT Infrared Deep Sky Survey (UKIDSS, Lawrence et al., 2007) Large Area Survey (LAS), which has a much smaller area than 2MASS but reaches limiting magnitudes of 20.5, 20.0, 18.8 and 18.4 mag in and bandpasses, respectively. The deepest near-IR photometry of Virgo galaxies corresponds to the Spectroscopic and H-band Imaging Survey of Virgo cluster galaxies (SHiVir, McDonald et al., 2011), that consists of pointed -band observations with the UH2.2 meter, UKIRT, and CFHT telescopes in good seeing conditions (). While SHiVir provides deep -band imaging for 286 VCC galaxies, it does not offer the one-to-one contiguous mapping with the optical NGVS data like NGVS-IR does.

One of the most immediate aims of NGVS-IR is to improve source classification in order to obtain the cleanest possible samples of objects of any given type, and of GCs in particular. For this purpose, it was decided to produce a first set of mosaic images and catalogs with a focus on compact sources with small angular sizes. As will be shown below, adding near-IR photometry to the optical data provides an impressive separation of Virgo GCs from foreground stars and background galaxies in near-UV/optical/near-IR color space. The diagram5, in particular, proves to be an exceptional tool for future searches of GCs and compact stellar systems, such as Ultra-Compact Dwarfs (UCDs, see Hilker, 2011, and references therein) around galaxies in a large variety of environments. With more than 20 deg of near-IR data soon available for Virgo, we will improve the decontamination of the optically-selected samples currently analyzed within the NGVS collaboration and provide extensive studies of the GC distributions in luminosity, color, and space. Those distributions contain a record of GC formation mechanisms and GC system assembly histories. They are also a key test for stellar population models, which still struggle to exploit optical and near-IR photometry jointly (e.g. Pessev et al., 2008; Taylor et al., 2011; Riffel et al., 2011). Finally, the construction of clean GC samples facilitates very efficient spectroscopic follow-up campaigns with little contamination by other sources. This, again, turns GCs into probes of galaxy cluster dynamics out to large scales (e.g. Schuberth et al., 2012; Romanowsky et al., 2012; Pota et al., 2013)

A further motivation for NGVS-IR is the study of UCDs and the dwarf galaxy – globular cluster transition. In the past decade, studies of the central regions of the Virgo and Fornax galaxy clusters have revealed dozens of massive compact objects in the magnitude range  mag with half-light radii between 10 and 30 pc (Haşegan et al., 2005; Jones et al., 2006; Misgeld & Hilker, 2011). Currently, two popular formation scenarios are that UCDs are the remnant nuclei of dwarf galaxies that were stripped by interactions (e.g. Bekki et al., 2003), or alternatively that they are agglomerates of multiple clusters that merged early in their lifetimes during violent star-formation episodes such as those triggered by galaxy mergers (e.g. Fellhauer & Kroupa, 2002). These scenarios can be tested through studies of the mass-metallicity relation (e.g. Harris et al., 2006). NGVS will provide us with a larger sample of massive compact stellar systems within Virgo than ever before, and with NGVS-IR we will be able to quantify the ages, metallicities and masses of these objects to higher accuracy than with optical colors alone (see e.g Puzia et al., 2002; Hempel & Kissler-Patig, 2004; Pessev et al., 2008).

The paper is organized as follows. We describe the details of our observations in section 2 and the data reduction steps in section 3. Sections 4 and 5 contain a discussion of the point source photometry and completeness estimates, while section 6 includes the presentation of our first results. We summarize our work in section 7.

## 2. Observations

### 2.1. Instrumental Setup

The NGVS-IR observations for the pilot field were carried out at the Canada-France-Hawaii Telescope (CFHT), using the Wide-field InfraRed Camera (WIRCam; Puget et al., 2004). WIRCam is mounted at the prime focus and consists of four cryogenically cooled HAWAII2-RG near-infrared (near-IR: ) HgCdTe detectors arranged in a mosaic, with typical interchip gaps of 45″ and a plate scale of 0.3″ per pixel (see Figure 1). The array is operated at Kelvin with negligible dark current (/sec) and readout noise. The field of view (FOV) of the full mosaic is about on the sky. The observations took place over a series of observing programs spanning the time period from December 2009 to July 2010 split over the French and Canadian TAC time allocation (see Table 1). They are arranged in a mosaic of WIRCam pointing or tiles, that cover degrees of the Virgo cluster around the central giant elliptical galaxy M 87 (see Figure 2), and match four pilot project CFHT/MegaCam pointings of the NGVS program (Ferrarese et al., 2012). All observations were carried out in Queue Service Observing mode (QSO) at an airmass less than 1.2 under clear observing conditions, when the estimated seeing was better than 0.8″.

### 2.2. The Dithering Pattern

#### General Considerations

Pointed near-IR observations are usually executed with a dither pattern that sequentially observes the science target followed by a blank-sky region to account for the strongly varying near-IR sky surface brightness. This strategy involves significant overhead and generates many sky observations that have generally no scientific use. In large surveys of non-crowded fields one may use the area of interest itself to estimate the sky, by implementing dithered observations. The dithering pattern also serves to fill the gaps between individual detector chips and defect areas of the detectors (of which WIRCam has quite a few). The NGVS-IR field, however, targets an area containing very extended galaxies. Surface brightness profile analyses of M87 found a (King, 1978; de Vaucouleurs et al., 1991; Liu et al., 2005), which is similar to the FOV of a single WIRCam chip. We developed a dedicated dithering strategy to obtain sky estimates within the survey area despite the presence of large galaxies. We use an individual integration time of 25 seconds to avoid saturation on the sky (the near-IR sky surface brightness on Mauna Kea10 varies around AB mag/arcsec). A combination of small and large offsets, the latter of which are defined based on archival images of the field, are used to sample the sky on a timescale shorter than the typical timescales of variations in its brightness.

#### The Exposure Dither

Each visit to one of the 36 pointings, or tiles, of the survey area is implemented as a sequence of 4 dithered 25 second exposures. This exposure dither is illustrated in Figure 2 by the WIRCam tile footprint that overlaps NGC 4486 (M87). The exposure dither pattern has a square geometry but this square is tilted with respect to edges of the detector. This facilitates the process of neutralizing smaller chip defects and covering the 45″ inter-chip gap, while keeping each sub-integration safely in the linear regime of the detector. It also keeps individual visits to a given pointing reasonably short compared to usual timescales of sky surface brightness variations ( in 10 minutes), so a given 4-exposure visit to one tile may be used to estimate the sky in tiles visited just before or just after.

#### Tile-to-Tile Sequencing

After one seconds visit is completed we always move the telescope to a new pointing in the 36-tile mosaic of the WIRcam survey. How these large offsets are chosen is based on the relative crowding of the tiles. Figure 3 shows the mosaic overplotted on a 2MASS archive image, together with a sample of photometrically selected GCs indicated as red dots (Peng et al. in preparation). Each tile is assigned a category A, B or C based on the crowding and spatial extent of objects that it covers. We use two types of sequences: the (A-C-A) observing sequence combines the most crowded tiles (C-type) with two non-crowded tiles (A-type), while the shorter (A-B) sequence combines one of the moderately-crowded tiles (B-type) with one non-crowded tile (A-type). The arrows in Figure 3 illustrate which crowded (B) and (C) tiles are observed together in sequence with (A) tiles which are relatively devoid of crowded regions and extended objects. This technique minimizes the required telescope slewing time and will allow accurate modeling of galaxy surface brightness profiles.

#### The Pointing Dither

Each (A-B) and (A-C-A) sequence is executed 27 times to reach the nominal exposure time of 45 minutes per pointing. We optimize the telescope nodding to simulate an on-tile dithering pattern. In other words, subsequent repeats of a given (A-B) or (A-C-A) are offset. In contrast to the exposure dither (see above), we call this pattern the pointing dither. It guarantees that each time all pixels cover a different section of the sky while the extent of the pointing dither allows us to combat the chip gaps and larger chip defects. The outer regions of the pointing dither pattern are covered by overlap regions between individual tiles of the entire WIRCam mosaic.

An optimal dither pattern needs to be compact, be scalable, and prevent redundant pointings. Such a pattern is realized in nature by the arrangement of flower leaves (phyllotaxis) and seeds on flower heads, a prominent example, for instance, being the sunflower head. Mathematically these patterns are described by a Fermat spiral (Vogel, 1979). This is a special type of the Archimedean spiral and can be analytically expressed as where and , while is the index of the individual pointing. The angle is the Golden Angle. Figure 4 shows the arrangement of the 27 pointings we use for moving between repeated (A-C-A) or (A-B) sequences. We tested several Archimedean spirals with different scalings and found that the Fermat spiral yields the least pixel-to-pixel variance in terms of x-y coordinate coverage (see the histograms in the left panel of Figure 4). In fact, the Fermat spiral pattern is the most compact and homogeneous coverage of a 2-D surface without x-y pointing redundancy. We scale it to the dimensions ( so that the central WIRCam inter-chip gap introduces the least variance on the final sky coverage map, which is illustrated in the right panel of Figure 4.

The pointing dither and the exposure dither ( seconds) work together to homogeneously fill the entire survey area and we refer, in the following, to the full dither sequence as F2D27. The individual components of F2D27 are illustrated in Figure 2 where the exposure dither is shown for tile #15 (center of the field), while the pointing dither is depicted for tile #1 (lower left corner of the mosaic). Each tile in the WIRCam mosaic is offset by 19′ from its neighbour in and to assure enough overlap between the individual mosaic tiles, so that the sampling of the WIRCam inter-chip gap remains the most critical for homogeneous sky modeling and similarly for the sampling of extended surface brightness structures, such as those of large galaxies. A coverage analysis of this WIRCam inter-chip gap region shows that our F2D27 pattern covers all sky locations more than 64 out of times, while the vast majority of the inter-chip gap is covered at least 80 times. This translates into a SNR fluctuation of only 14% and is fully propagated by our pipeline to the variance maps.

In summary, the observing strategy includes the following exposure time outline for the two types of sequences:

 \centering(A1−C−A2)\@add@centering → 3⋅27×(4×25)sec=8100sec (A−B) → 2⋅27×(4×25)sec=5400sec

and contains a total of six (A-C-A) sequences and nine (A-B) sequences that cover the entire field. Due to observing scheduling constraints, sequence (A6-B5) could not be observed and thus will not be discussed in this work. On the other hand, some of the observed sequences benefited from more than the requested 27 visits.

## 3. Data Reduction

The data reduction process was divided into three major parts: the pre-processing stage, the main image processing, and the image post-processing.

### 3.1. Pre-Processing

The processing of the raw data was done using the iwi pipeline v2.011. The pipeline involves the following steps: flagging saturated pixels, non-linearity correction, dark subtraction, flat fielding and bad-pixel masking.

The saturated pixels were identified in the raw images using a threshold value of 36000 ADU. A record of these pixels was kept and they were flagged with the value 65535 in the final pre-processed image. The non-linearity of the WIRCAM detectors is about 5% at 30000 ADU and it was corrected at the beginning of the pre-processing. Dark frames were obtained at the beginning and the end of each observing night, and a master dark frame was computed by taking the median of 15 dark frames. The dark subtraction has a negligible impact on the pre-processed images, since the dark current of WIRCam detectors is significantly less than /sec.

Twilight flat fielding was applied by the iwi pipeline: Twilight flats are obtained at the beginning of each observing night and a daily twilight flat is built by taking the median of 15 consecutive frames. The pipeline divides the semester into several sets of observing nights and then computes a master twilight flat on each set. All science images were corrected by dividing by the corresponding normalized master flat.

The master bad pixel masks were built by analyzing the normalized master flats previously computed. Bad pixels were identified by applying a sigma threshold algorithm. The bad pixel masks were constructed to have a value of 1 for good pixels and 0 for bad pixels.

We refer to the pre-processed images in the following as science images.

### 3.2. Main Image Processing

The processing consisted of masking cosmic-rays and satellite-trails in the science images, removing the residuals left by saturated stars, subtracting the sky from the images, computing the astrometric solution, applying the photometric calibration and producing the final stacked images.

#### Cosmic-Ray and Satellite-Trail Removal

The cosmic rays were identified by running the Laplacian cosmic ray identification algorithm (LACosmic; van Dokkum, 2001) on the science images. This algorithm is based on a variation of the Laplacian edge detection method and cosmic rays are identified by the sharpness of their edges. The code convolves the image with a laplacian kernel, then identifies cosmic rays as sharp borders in the image and creates a cosmic-ray mask for each science image. The best LACosmic parameters were obtained by performing a visual inspection of the masks and verifying that the peaks of non-saturated bright stars were not miss-identified as cosmic rays.

The detection of satellite tracks is not possible in the pre-processed images, since the sky brightness dominates the image while a typical satellite track is hundred times fainter than the sky. We run a quick sky removal on each image (target) that consisted of identifying the 10 closest images in observing time (sky), then taking the median of the sky images pixel by pixel and subtracting it from the target image. The sky subtracted images were smoothed with a boxcar average of width 20 pixels and the coordinates of pixels with signal-to-noise ratio greater than 1.2 were registered. Satellite tracks were identified as straight lines in the XY detector plane of the registered coordinates.

#### Saturation Correction

Bright and saturated stars leave residuals in the science images that last for about six subsequent images. First, the saturated bright stars were identified as those regions of at least four connected pixels with counts higher than 65535 and not masked as bad pixels in the master bad pixel mask. For chips 1 and 3, it was necessary to mask the subsequent three images, while for chips 2 and 4 it was necessary to mask the subsequent six images due to the different chip and detector characteristics. The science images were sorted according to their observing time, the pixels belonging to a saturated star were identified and the six (or three) consecutive images after saturation were masked at the corresponding coordinates.

#### Sky Subtraction

One motivation for the (A-C-A) or (A-B) observational strategy presented in Section 2.2 was to design a high-quality sky subtraction method. The observations in each (A-C-A) or (A-B) sequence were declared as target or sky tiles, and for each target image there was a set of sky images defined for computing and subtracting the sky in the corresponding target image. The sky images from A-type tiles are generally less crowded regions than target images, but do still contain a significant number of point-like and extended sources that need to be masked before computing the final sky image. The sky subtraction method is a two-step iteration process: the first iteration consists of computing a median sky image for each target image, then subtracting it from the target image and finally building a stacked target image from all such treated frames in a sequence; the second iteration consists of identifying the sources in the stacked target image, masking these sources in the each individual sky image, computing a new median sky and finally subtracting it from the target image.

The selection of the set of sky images to model and subtract the sky contribution in each target image depends on the amplitude and variance of the sky surface brightness and is based on the following three criteria:

the category of the (A-C-A) or (A-B) observing sequence that the target image is part of,

the time difference between sky images and the target image time stamp,

the need for a sufficient number of sky images to compute a reliable median sky image.

The first criterion is related to the sky-target dithering sequence explained in Section 2.2.3, and consists of the following set of rules: 1) if the science image belongs to type-A tile, then sky images are selected from type-A tiles of the same sequence; 2) if the science image belongs to a type-B tile, then sky images can be selected from type-A and type-B tiles; 3) if the science image belongs to a type-C tile, then sky images can be selected only from type-A tiles.

The second criterion consists of choosing sky images close enough in time to the target image. We used only those sky images taken within equal or less than 15 minutes before and after the corresponding target image observation. In particular, at the beginning and the end of each observing run there were sets of sky images containing just 8 images, and as it will be explained by the third criterion, that number of images is below the recommend number of frames to model the sky.

The third criterion consists of defining the sufficient number of sky images that is necessary to compute a reliable median sky image that robustly captures the sky variations. We ran several quality tests on the sky images and the conclusion was that the minimum number of frames for this survey is ten sky images given the variability of the near-IR sky throughout our observing campaign. This criterion was applied as follows: we sorted the sky images according to their time stamp with respect to the target image, then we computed the maximum between 10 and the number of sky frames that fulfilled the first and second criterion, and the final set of sky images was defined by applying that maximum number.

Once the set of sky images is defined for each target image, we compute the pixel-to-pixel median sky image and subtract it from the corresponding target image.

### 3.3. Post-Processing

We applied post-processing routines to the images because some of the processed images were still showing visible systematics. These post-processing steps consist of removing the slope from each amplifier detector and subtracting large-scale variations of the sky background.

#### Residual Amplifier Differences

Each WIRCam detector has 32 amplifiers with different gains that result in background level offsets that differ by in the twilight flat-field images. Most of the signal of the amplifiers disappears after flat-field correction and removing the sky during the pre-processing pass, but there is still a non-negligible contribution that appears as horizontal bands through the image. To model and subtract this residual, we start with sky-subtracted images obtained from the processing pass and select the pixels belonging to each amplifier, taking the median along the x-axis of the detector. Then we fit the variation with a linear relation as a function of the y-axis coordinate and finally subtract this model from each horizontal band in each of the sky-subtracted images.

#### Large-Scale Variations

After close examination of the gain corrected and sky subtracted images, we found large-scale variations of the image background still present after the sky removal, with an amplitude of 0.3% of the typical sky. These feature are likely to be produced by a variation in the thermal stability of the detector and/or internal reflections within the instrument due to the presence of large galaxies. They have a non-negligible effect on the stacked the images.

In order to improve the background flatness, we mask the sources and model the large-scale background structure in each frame. For this purpose, we first use Swarp (Bertin et al., 2002) to build a stacked image using all the sky-subtracted frames computed in the processing stage. We run SExtractor on the stacked images (v2.5.0, Bertin & Arnouts, 1996) and obtain the segmentation maps with pixels associated with the sources. To remove the outer low-surface brightness parts of galaxies, we increase the sizes of the masked regions by 50 to 800 pixels, depending on the size of the galaxy, and project the source mask back into the individual pre-stack sky-subtracted images. The large-scale variations are computed on each of the four WIRCam detectors independently. Each detector is divided into a grid of cells. Those cells in which more than 50% of the pixels are masked are discarded. In all other cells a median background value is computed. A large-scale variation image is then computed using the Kriging linear interpolation algorithm as implemented in IDL (v8). This background model is then subtracted to produce the final science-grade images. This method has some limitations when dealing with galaxies with extremely extended surface brightness profiles. However, these limitations have no impact on the analyses and results presented in this paper.

### 3.4. Astrometric Calibration

We calculate the relative astrometric solution with the Scamp software package (v1.7, Bertin, 2006), which reads the catalogs produced by SExtractor and cross-matches the source positions against an astrometric reference catalog, such as 2MASS (Skrutskie et al., 2006). We run Scamp using a maximum search range of ′ (position_maxerr=0.15) and a minimum signal-to-noise ratio of 40 (sn_thresholds = 40, 80) to match reference stars.

We compute the distortion maps for each of the WIRCam detector chips and find that the pixel scale varies in a radially symmetric way by at most 0.5% between the center of the detector mosaic and the outer radius of the WIRCam FOV. The astrometric calibration against the 2MASS reference frame is based on  stars and the resulting astrometric world coordinate solution accuracy has an approximately gaussian distribution with a FWHM smaller than 0.02″.

### 3.5. Photometric Calibration

We use the 2MASS point source catalog (2MASS-PSC; Cutri et al., 2003) to perform the preliminary photometric calibration of the sky subtracted images. We select 2MASS point-like sources with photometric quality flags equal to A and B (i.e. reliable photometry) that correspond to sources with S/N ratio greater than 7 and magnitude error less than 0.16 mag.

The photometric calibration was done using Scamp as for the astrometric calibration. Scamp distinguishes between two types of exposures: those observed under photometric conditions and the others. The information about the photometric nights was obtained from the QSO weather reports and the images observed during those nights were visually inspected to define a clean sample of exposures taken in photometric conditions. The zero-point of the photometric calibration was computed using

 mstd = minst+zp−kX, (1)

where is the 2MASS magnitude, is the instrumental magnitude measured on the sky-subtracted images, is the zero-point of the photometric calibration, is the airmass term, and is the airmass of the observation.

The four WIRCam detectors have different gains and in principle a should be computed per exposure and detector, but the number of reliable 2MASS stars per detector can be as low as ten. We opt for computing the offsets between the WIRCam detectors using groups of exposures, and then keep these offsets fixed to assess variations between exposures in the group. The offsets are computed by grouping the exposures every 15 minutes, then crossmatching the detected sources with reliable 2MASS point-like sources and computing the median value of for all the stars on each detector. Figure 5 shows the distribution of the difference between the values of chip 1 and and those of chips 2, 3 and 4 ( offset), where the distribution for chip 2 is found to be the most concentrated. The offsets for detectors 2, 3 and 4 are 0.12, 0.08 and 0.03 mag, respectively.

For computing the airmass term () in Equation 1, we first compute the mean value of on each exposure and then apply a linear regression analysis using only photometric nights. The estimated value of for this survey is . We use Equation 1 to estimate the for exposures observed in photometric nights and add the value as phot_zp to the header.

For exposures observed in non-photometric nights, Scamp runs an internal crossmatch of point-like sources and then adjusts the for those exposures by minimizing Equation 2. We adopt the same parameters as in the previous section, but this time we use a S/N ratio of 20 or higher in order to have a larger sample of stars for doing the photometric calibration (sn_thresholds=20). We compute

 χ2=∑s∑a∑b>awab(zpa+minst,a−zpb+minst,b)2 (2)

where and refer to exposures that contain the point-like source in the overlapping region, and is the non-zero weight for the pair of detections in exposures and , which are computed via where is the square sum of all contributing photometric variances.

The instrumental magnitude, , in Equation 1 is measured on the sky-subtracted images with a circular aperture photometry of radius , where is the mean value of the full-width-at-half-maximum of the bright stars in each exposure. Figure 6 shows the distribution of the FWHM of point-like sources of all the images of the survey. The FWHM distribution is not homogeneous and covers the range between 0.4″ and 1.0″. Thus, the photometry aperture is large enough to include more than 95% of the total light. The photometric zero points adopted in this work are shown in Table 2. The refers to the mean zero point of science frames observed in the respective night and the quoted uncertainty to its standard deviation.

### 3.6. Stacking and Data Quality Assessment

The sky-subtracted images were stacked using the Swarp software (v2.19, Bertin et al., 2002) after the astrometric and photometric calibrations were completed. Swarp uses the astrometric solution calculated by Scamp and combines several frames to produce a final image of a given size and pixel scale. We use the -band quadrants of the NGVS pilot field area as a reference for building the same size field-of-view -band stacked images, since it helps the subsequent analysis. We note that the pixel scale of WIRCam images is 0.3″, while the pixel size of our final stacked images is 0.186″ and that we build the stacked images using only images with FWHM smaller than 0.7″ (see Figure 6). The median seeing of all selected data is 0.54″. The NGVS-IR pilot field covers a total area of 3.98 degrees, and consists of 34 slightly overlapping tiles with coverage of about 26.5′26.5′. For data handling purposes, we sub-divide the pilot field area into four quadrants roughly corresponding to about 1 degrees each.

The large majority of the sky subtracted images in the survey are super-critically sampled (above the Nyquist limit). We use the Lanczos-2 resampling algorithm which offers a good compromise between improving the image resolution and minimizing the number of artifacts introduced by the resampling routine. The images are combined using a sigma clipping algorithm that computes the dispersion on a pixel-by-pixel basis and rejects pixels beyond of the mean value.

The NGVS-IR image stacks, which this paper is based on, reach -band surface brightnesses fainter than AB mag/arcsec. In some of our NGVS-IR stacks we can make out dwarf galaxies with AB mag/arcsec. The extended emission of 94 % of the known VCC galaxies can be visually seen in the final images. However, we emphasize that the data products discussed here were not obtained with the purpose of optimizing surface brightness sensitivity and we, therefore, postpone any study of luminosity profiles to future work. Specific data reduction procedures will also be required for the study of near-IR surface brightness fluctuations in Virgo galaxies. In Figure 7 we illustrate the data quality with a color stack of an cutout around the interacting galaxy pair NGC 4435+4438, which includes five of the faintest VCC dwarf galaxies and even fainter low-surface brightness dwarfs all of which are detected on the stacks.

Before performing detailed completeness test, we note that visual inspection of the images shows that our NGVS-IR WIRCam survey achieves a point-source sensitivity of  mag, thus reaching about 2.9 mag deeper in this passband than the UKIDSS Large Area Survey, the most extensive and deepest recent survey covering the NGVS pilot field (Lawrence et al., 2007, 2012). NGVS-IR compares in surface brightness sensitivity with the pointed -band observations of the SHiVir survey12 of bright Virgo galaxies of (McDonald et al., 2011). These, however, are not contiguous images and were taken with several instruments in varying observing conditions.

## 4. Photometry

Source detection in the stacks and photometric measurements were performed with SExtractor (v2.5.0 Bertin & Arnouts, 1996) and PSFex (Bertin, 2011). After a first SExtractor run for source detection, PSFex produces local models of the point spread function (PSF), that can then be used in a second run of SExtractor to obtain integrated PSF magnitudes of stars. As input for PSFex we use a clean sample of stellar sources, obtained by combining a compactness and a color criterion (see Sect. 6.3), and by rejecting any star that saturated the detector in one or more individual images that entered the stack. The PSF was modeled separately across the four individual quadrants of the pilot field (see Figure 8). Each of these quadrants consists of WIRCam pointings, and PSF variations are dominated by pointing-to-pointing variations in the corresponding FWHM distribution functions (see Figure 6). The spatial variations of the PSF are modeled with polynomials of degree seven to allow for changes on the typical scale resulting from this pattern.

With the SExtractor parameters detect_thresh=1.4, analysis_thresh=2.0, detect_minarea=6 and deblend_mincont=0.000001, a total of 450,000 sources are detected and about 17,800 of these are identified as stars. The accuracy of our measurements are best described by focusing on these point sources. In Figure 9 we compare the WIRCam photometry with the corresponding 2MASS and UKIDSS data. The dispersion in both panels is predominantly due to the photometric errors in 2MASS and UKIDSS, respectively. The UKIDSS data (DR8, Lawrence et al., 2007) provides a tighter calibration relation than 2MASS. Both panels exhibit a significant color term. The left panel of Figure 10 shows that the color term relating WIRCam and UKIDSS magnitudes seems to depend on the -band magnitude. At the bright end, we see the actual effect of the different shapes of the filters of WIRCam and UKIDSS, while at the faint end the slope mostly reflects large uncertainties in the UKIDSS photometry. For the UKIDSS native magnitudes (in the Vega photometric system, see Appendix A.1) we have:

 (K−KUKIDSS)Vega=−0.09+0.27⋅(H−K)UKIDSSforK<17.0mag

The right panel of Figure 10 illustrates the external photometric accuracy after this simple color-term correction. The photometric comparison between UKIDSS and WIRCam data is affected by a statistical uncertainty of 0.05 mag down to about Vega mag, which corresponds to AB mag.

## 5. Completeness Estimates

The photometric completeness of point-sources in the final image stacks was computed by adding artificial stars with a range of magnitudes at randomly selected positions into the image. We use a set of IDL and Python scripts to generate a list with the positions and magnitudes of artificial stars, avoiding areas with the presence of saturated stars and very extended galaxies, e.g. the central regions of M87.

The first step consists of generating the positions of artificial stars. We run SExtractor on the final images and generate the segmentation maps. We use these maps to avoid adding stars on top of a galaxy or other high density regions. Then, the image is divided into a 2D grid of 200-pixel wide bins and we randomly generate 10 positions in each bin. For each one of these positions we check that there is no collision with other sources (real or artificial) in a radius of 16 pixels, thereby avoiding artificial crowding effects. The final list consists of about 100,000 sources for each final quadrant image.

In the second step we generate the magnitudes of the artificial stars. First, we run PSFex on the final images and compute a position-dependent PSF model. We use this PSF model to properly add a point-like source in different positions on the image. Then, we generate a homogeneous distribution of magnitudes in the range mag and assign them to the stars generated in the previous step ().

The third step adds the artificial stars to the image at the corresponding positions, with the generated magnitudes using the appropriate PSF model. We use a highly efficient multi-thread Python script for reading the list of coordinates and building the point-like sources according to the PSF model computed with PSFex using all available CPU cores.

We run 200 realizations of this procedure for each NGVS-IR quadrant ending up with artificial objects per quadrant. We run SExtractor with the same parameters used for the original stacks (see Section 4) and cross-match the positions of the detected sources with the artificial stars generated in step one, hereafter called recovered stars. The top panel of Figure 11 shows the difference between the magnitude of artificial stars in tile +0+0 as measured with SExtractor () and assigned in the mock catalog (). This figures illustrates the internal photometric errors of the NGVS-IR survey and demonstrate the excellent quality of the stacked images. In order to estimate the completeness, we adopt a magnitude bin of 0.1 mag and compute the ratio between the number of recovered stars and the number of artificial stars as a function of magnitude and position. The bottom panel of Figure 11 shows the completeness fraction as a function of the input AB magnitude for each of the four quadrants, the entire survey, and the deepest and shallowest regions in the pilot field. The 90% completeness magnitude for point-source detections and the corresponding limiting magnitude together with a conservative surface brightness limit estimate for the NGVS-IR quadrants are summarized in Table 3.

Figure 12 illustrates the variations of all the completeness magnitudes across the four quadrants. The main reason for the fluctuations visible across the field stems from the photometric quality, i.e. FWHM distribution (see Figure 6) in each individual WIRCam tile contributing to the final quadrant stack.

## 6. Discussion

### 6.1. Color-Color Diagrams

Adding more filters and a successively wider SED coverage to the investigation of astronomical objects adds more diagnostic power to the analysis and delivers ultimately more robust and astrophysically meaningful results (e.g. Park & Choi, 2005; Puzia et al., 2007). Color-color planes are efficient tools for the classification of sources in large-scale imaging surveys (e.g. Daddi et al., 2004; Faber et al., 2007). By combining near-UV and optical data from NGVS and near-IR photometry from the NGVS-IR, the advantage of spanning the entire spectral range of stellar emission from the atmospheric UV cutoff at  Å to the near-infrared at m becomes clear when one is confronted with the vastly improved system throughputs and SED coverage of the combined NGVS+NGVS-IR filter set as illustrated in Figure 13.

To construct a first NGVS+NGVS-IR source catalog, we cross-matched the NGVS-IR catalog positions with those of the NGVS catalog, which contains approximately half a million objects. The distribution of angular separations between and -band coordinates peaks at 0.09″, and 95% of the detected sources (of small angular size) have a separation smaller than 0.3″. This small bias is due to the marginally resolved nature of background galaxies which constitute a large fraction of the matched catalog and introduce a random offset in the matching accuracy that appears as a systematic in one radial dimension. Future dedicated versions of those matched catalogs will be built with object-dependent procedures that allow for self-consistent optical and near-IR measurements, in particular for extended sources. However, independent measurements and subsequent matching are fully sufficient for the qualitative purposes of this paper.

Historically, the optical/near-IR color-color plane that is most widely used in conjunction with large extragalactic surveys is the diagram which combines and color indices or equivalents thereof. The diagram was highlighted by Daddi et al. (2004) as a means of identifying both passively evolving and star-forming galaxies located at redshifts larger than 1.4 or so. It was published later for a variety of deep fields, confirming the original segregation between low and high-redshift sources (e.g. Lane et al., 2007; McCracken et al., 2010, 2012; Bielby et al., 2012). The technique is now being frequently used for studies of the clustering properties, intrinsic morphologies, stellar populations, and the X-ray luminosities of -selected galaxies (e.g. Lin et al., 2012; Yuma et al., 2011, 2012; Ly et al., 2012; Rangel et al., 2013). Ultra-deep follow-up spectroscopy campaigns confirm the efficient selection of this type of high-redshift, star-forming galaxies that suffers contamination by other sources (Kurk et al., 2013).

### 6.2. The NGVS+NGVS-IR analog of BzK: the gzKs color-color plane

In terms of SED coverage of the involved filters, the closest analog to the diagram that can be constructed from our combined NGVS+NGVS-IR data is the diagram which is shown in the top panel of Figure 14, which illustrates our pilot field sample greyscale-coded by the Gaussian propagated total error of each contributing filter (see also Figure 13 for a comparison of the corresponding vs.  system throughput curves). The mean foreground extinction towards M87 is , , and mag, and was taken from Schlafly & Finkbeiner (2011)13. The following analysis includes the corresponding extinction correction terms. We observe various characteristic object overdensities and sequences in the plane, such as the narrow sequence of foreground Galactic stars and other less constrained object distributions including actively star-forming and passively evolving galaxies at various redshifts that were discussed in previous works (e.g. Bielby et al., 2012; Merson et al., 2013). We use equations (1)–(8) provided in Bielby et al. (2012) to transform the object selections into our color-color plane. In particular, we use the median color mag of our entire sample and obtain for the selection of star-forming galaxies at the following relation (shown as solid magenta line in Figure 14):

 (z−Ks)0>1.233×(g−z)0−0.017. (3)

To select passively evolving galaxies at we obtain (illustrated as dashed magenta line in Figure 14):

 (z−Ks)0<1.233×(g−z)0−0.017∩(z−Ks)0>2.5. (4)

To select foreground stars and separate them from galaxies our transformations yield (dotted magenta line in Figure 14):

 (z−Ks)0<0.37×(g−z)0+0.474 (5)

The corresponding relations shown in the top panel of Figure 14 illustrate the quality of separating the classically defined galaxies at from the general locus of more nearby galaxies and the stellar sequence. However, it is quite evident that even with the NGVS+NGVS-IR high-quality data, a clear separation of stellar foreground from the background galaxy population is not entirely possible based on the color-color diagram alone.

Overall the color-color plane is the tool of choice for the study of the redshifted universe beyond Virgo due to the superior photometric depth of the filter compared with the band observations. However, unlike any other deep optical/near-IR survey, NGVS+NGVS-IR also contains a nearby structure: Virgo itself. The most striking difference between our color-color diagrams and those of other surveys is not due to Virgo galaxies – these occur in negligible numbers in any given part of the diagram –, but due to GCs that are concentrated around the massive elliptical galaxy NGC 4486 (M87) and other giant ellipticals located in the central regions of Virgo (see Figure 3). To demonstrate the locus of colors of Virgo GCs we highlight radial-velocity confirmed GCs in Figure 14 as red dots (E. Peng et al., in preparation). This GC locus overlaps with the colors of stars at the blue end, and is heavily contaminated by background galaxies towards redder colors, in the diagram.

### 6.3. The uiKs Color-Color Plane

For a photometry-based selection of GCs and other compact stellar systems, such as UCDs and dwarf galaxies, we show in the following that the most powerful color combination is given by the diagram, which combines high-quality near-UV, optical, and near-IR photometry from NGVS and NGVS-IR. The bottom panel of Figure 14 shows the color-color diagram, in which the data is deredden with the values extracted from the Schlafly & Finkbeiner (2011) maps using , , and mag. The typical structures seen in the diagram (top panel of Figure 14) that were classified in previous deep-field surveys appear much more prominent and better defined in the color-color plane. The stellar sequence seen in the diagram remains very clearly identified in this new plane and appears more separated from the main cloud of galaxies. This will be particularly useful for the analysis of stellar age and metallicity distribution functions in the Virgo overdensity described in Ferrarese et al. (2012). Several new narrowly defined features in the plane become visible at intermediate colors that we attempt to classify qualitatively in the following. Most importantly, however, we note that the contamination of the GC locus is very significantly reduced with the use of the filter combination.

#### An Efficient Tool for Star Cluster Selection

In order to better understand the various features in the plane, we overplot in Figure 15 predictions of simple stellar population (SSP) model calculations based on a customized version of the population synthesis code Pégase (Fioc & Rocca-Volmerange, 1997) that includes the exactly matched throughput functions for all NGVS+NGVS-IR filters. While in the plot these SSP models coincide with the overlap region between stars and galaxies, in the plane they fall right on top of a sharply defined sequence which we identify as GCs. We note that the metallicity and age coverage of the shown SSP model predictions, i.e.  and to 13 Gyr (each running from bluer to redder colors), agrees fairly well with the GC candidate locus, in particular in the plane, and is consistent with what is expected from previous photometric and spectroscopic studies of stellar populations in typical Virgo GCs located in the vicinity of M87 (e.g. Hanes & Brodie, 1986; Cohen et al., 1998; Hanes et al., 2001; Kissler-Patig et al., 2002; Jordán et al., 2002; Tamura et al., 2006b; Peng et al., 2009; Yoon et al., 2011; Forte et al., 2013). We verify the fidelity of the GC selection with objects that have the systemic radial velocity of the Virgo cluster using spectroscopy available in the literature (Hanes et al., 2001; Strader et al., 2011) and recently obtained with the multi-object, moderate-dispersion spectrograph Hectospec at MMT (E. Peng et al., in preparation). In Figure 14 we overplot such radial-velocity confirmed Virgo GCs in the and color-color diagrams as red symbols. At this point we defer a more quantitative analysis of the Virgo GC stellar populations to a future paper and continue with the investigation of the other sequences, keeping in mind the -based photometric GC selection.

#### The Redshift Evolution of Galaxies in the uiKs Plane

In addition to the SSP model color predictions, we plot in Figure 15 the redshift evolution in the and color-color space of three prototypical SEDs of galaxies born at with different types of star-formation histories. These tracks were computed using the NGVS+NGVS-IR instrument throughput curves (see Figure 13) with a customized version of the population synthesis code Pégase (Fioc & Rocca-Volmerange, 1997), but other population synthesis codes show very similar trends (see also Dahlen et al., 2013). In the order of numerically increasing colors, the tracks correspond to 1) a continuously star-forming galaxy, 2) a galaxy that formed the majority of its stellar content at high redshift and star formation continues at a lower constant rate, 3) a passively evolving galaxy that formed all its stellar content at .

We notice that, in contrast to the plane, galaxies with any level of star formation are expected to be fairly distinct in their color properties from GCs and foreground stars from out to redshifts and higher. We find that the hook in the redshift evolution of star-bursting galaxies at in the color-color plane is responsible for much of the contamination at colors and mag, which is the mean locus of metal-poor and intermediate-metallicity GCs. Such a contamination is entirely avoided in the plane which is due to the mostly blue-ward evolution of the color when redshifting the steeply increasing near-UV part of star-forming galaxy SEDs. The redshift evolution of the corresponding passively evolving SEDs overlaps with the local galaxy background of Virgo and its diagnostic power is only limited by the depth of our -band photometry. The shallower magnitude limit of the -band data is apparent when the blue edge of the galaxy distribution in each diagram of Figure 15 is compared. Although the redshift distribution of the background galaxies in the plane cuts off at lower values as compared to the plane, especially along the sequence of passively evolving galaxies (the cut-off is near instead of near ), the plane provides overall a clearly linear and powerful selection diagnostic of star-forming galaxies in the redshift range . Because of this property and the photometric depth of the NGVS+NGVS-IR data we are in the position of searching for galaxy clusters at redshifts , which will be described in detail in future papers.

#### The Diagnostic Power of the uiKs Plane

The reason for the superior GC-star and GC-galaxy separation power of the plane lies in the combination of the near-UV, optical, and near-IR color information, compared to the classic optical/near-IR plot (see Fig. 14). While the optical/near-IR color probes the effective temperature of the red giant branch (RGB) stellar population, the near-UV/optical color probes the stellar fluxes blue- and redward of the  Å break and is, therefore, sensitive to hot stellar evolutionary phases. This leads to a very efficient GC-galaxy separation due to the sensitivity of the -band to hot stellar components in star-forming galaxies, combined with the redshift evolution of galaxy SEDs (see Figs. 13 and 15). The efficient GC-star separation in the diagram is due to the fact that at a given color (i.e. effective RGB temperature), the -band picks up the additional hot stellar component flux generated by horizontal branch stars, making GCs appear bluer in the color.

It is worth pointing out that attempts at photometrically separating extragalactic GCs from foreground stars and background galaxies exist in the recent literature, however, based on a combination of purely optical/near-IR (Puzia et al., 2002) or purely near-UV/optical colors (Kim et al., 2013). Neither of these color-color planes can separate foreground stars, GCs, and background galaxies. The only earlier study that combined near-UV/optical/near-IR photometry of relatively small GC samples was attempting to study the age distribution function of GCs in giant elliptical galaxies (Hempel & Kissler-Patig, 2004), and was based on much poorer photometric quality. To our knowledge, the diagnostic plane presented in Figure 14 is the first such plot based on high-quality near-UV/optical/near-IR photometry that allows for an efficient GC-star and GC-galaxy separation.

In order to assess the required photometric quality for robust GC-star and GC-galaxy separation, we plot in Figure 16 the and color-color plane and parametrize the symbol colors by the photometric uncertainty in each contributing filter. It becomes immediately clear that the photometric depth, i.e. photometric error distribution function in each filter, influences in a non-linear, and sometimes quite dramatic way the selection of galaxies as defined in Daddi et al. (2004). The corresponding relations are shown in the left panels (as in Figure 14) and illustrate the quality of selecting the classically defined galaxy samples as a function of photometric depth. It is also quite evident that even with the NGVS+NGVS-IR high-quality data, a clear separation of stellar foreground from the background galaxy population is not entirely possible based on the color-color diagram alone, in particular Virgo GCs cannot be robustly identified.

The right panels of Figure 16 demonstrate the high potential for studies of stellar populations in the nearby universe. We find that datasets with a photometric accuracy of , , and mag (corresponding to , , AB mag in our dataset) will have the potential to select very clean samples of extragalactic GCs. We will discuss and quantify the exact contamination fractions of such samples in future papers.

With the superior spatial resolution of our -band data (median seeing , see Section 3.6 and Figure 6) we search for additional diagnostic tools related to the morphology of objects in our survey field. For this, we devise a simple magnitude difference between a 3″ diameter aperture magnitude and a PSF magnitude. The PSF magnitude produces integrated colors for stars very similar to the total integrated colors for GCs, and a color roughly representative of the central colors for galaxies. The results for all our objects with -band photometry are shown in Figure 17 as a color parametrization of the plane according to the morphological measure of their radially symmetric compactness. Intriguing correlations between this simple morphology parameter and colors appear, that will be studied in subsequent papers of this series. We note here, that the figure clearly shows that GCs are marginally resolved in the WIRCam images: their observed shapes are slightly less concentrated than those of stars. The accuracy of measuring relative GC sizes is sufficient to even make out the classic size difference between blue and red GCs (e.g. Webb et al., 2012, and references therein). Attempts to identify GCs solely on the basis of their morphology, however, produce samples whose colors reveal a high level of contamination by stars and remote galaxies. Our final classification algorithm (in preparation) will therefore be based on both colors and morphology.

## 7. Summary

In this paper, we have described the NGVS-IR survey, its motivation, observational strategy, and science goals. We discuss in detail the data reduction procedures and present -band imaging data of the central 4 deg of the Virgo galaxy cluster observed with WIRCam mounted on the Canada-France-Hawaii Telescope (CFHT). These NGVS-IR images have a median seeing of and a 90% completeness magnitude of better than AB mag, thus, superseding in all aspects the 2MASS and UKIDSS mosaics of the same region.

Using the full near-UV to near-IR SED coverage of our NGVS+NGVS-IR data we investigate the color-color plane as a diagnostic tool to identify and select characteristic object classes. We also study the color-color plane as the closest equivalent to the plot from which star-forming galaxies at redshift are typically selected. With the described diagnostics we can identify several distinct and isolated groups of objects in the diagram that correspond to i) star-forming galaxies out to redshifts and higher, ii) a mixed distribution of galaxies with a variety of (non-negligible) levels of on-going star formation, iii) a redshift sequence of passively evolving old galaxies, iv) a relatively smooth sequence of Virgo GCs and UCDs, and v) a very distinct sequence of foreground Milky Way stars.

According to the photometric error distributions of our NGVS+NGVS-IR data in the color-color diagram (see right panels of Figure 16), we can isolate the Virgo GC population, virtually free of contamination from foreground stars and background galaxies, if we restrict our selection to objects with photometric errors , and mag. Independent of the photometric quality of the dataset such a selection would be fundamentally impossible from the diagram. The diagram is to our knowledge the most efficient and least biased extragalactic GC selection technique that is relying purely on photometric colors and is likely to prove extremely efficient in the study of distant star cluster systems with the combined near-UV/optical/near-IR instrumentation that will be available on the upcoming generation of survey telescopes, such as the Large Synoptic Survey Telescope (LSST) and the Euclid spacecraft, as well as future facilities, like for instance E-ELT, JWST, etc.

## Appendix A WIRCam Ks photometry - additional information

### a.1. AB magnitudes versus Vega-based magnitudes

The constant difference between Vega-based magnitudes and AB magnitudes for the WIRCam Ks observations can be obtained from synthetic photometry. The CFHT WIRCam web pages14 warn that various authors have used conversion constants that differ by up to 0.2 magnitudes. How the various values have been obtained, however, is not documented systematically and is rarely indicated in journal publications. We detail our own derivation below.

Figure 18 shows the components of the WIRCam Ks transmission curve. The data for the instrument components were obtained from the CFHT WIRCam Throughput web pages (files cfh8302.dat for the filter, and file WIRCamOpticsResponseCurve.xls for the optics). The optics transmission curve includes all optics except the telescope mirror and tip-tilt plate. We thus implicitely assume these two components have a flat response. The atmospheric transmission for Mauna Kea was obtained from the website of the Gemini Observatory. It is based on the ATRAN code (Lord, s.D. 1992, NASA Technical Memorandum 103957) and assumes a water vapor column of 1mm.

Table 4 provides the conversion constants obtained for WIRCam Ks under various assumptions. Three models for Vega were used. alpha.lyr.stis.003 and alpha.lyr.stis.005 are two versions of the standard Vega spectrum distributed by the Space Telescope Science Institute (Bohlin 2007 ASP Conf Ser. 364, 315, and calspec webpage on www.stsci.edu). They differ only shortward of 0.53 m, which affects the NGVS and bands but not WIRCam Ks. VegaLCB is the reference spectrum distributed with the population synthesis code Pégase (Fioc & Rocca-Volmerange 1997). The differences with the HST spectra are (1) a globally lower flux level (0.7 %), (2) a lower spectral resolution (which leads to small differences of the integrals over hydrogen features), and (3) local differences in the shape of the spectrum, such as between 2.3 and 2.5 m or around 1.5 m. The conversion constant we recommend is

 Ks(AB)−Ks(Vega)=1.827mag

This value is essentially identical to the value obtained by Stéphane Arnouts at CFHT (1.824 mag).

Users of Pégase should be reminded that code provides Vega-type magnitudes and colors under the assumption that Vega has a magnitude of 0.03 (and not 0.0) in all filters. Therefore, the WIRCam Ks magnitudes produced for synthetic populations by Pégase in the AB-system and in the Vega-based system, using VegaLCB as a reference, differ by 1.804 instead of 1.834. Rather than using Vega, our customized version of Pégase computes AB synthetic magnitudes directly.

## Appendix B Extinction coefficients

We have derived extinction coefficients by comparing the synthetic colors of stellar model spectra and of reddened versions thereof, using the extinction law of Cardelli et al. (1989) with R=3.1. The stellar models used were representations of Vega, of the Sun and of a red star (in practice, the solar model seen through 3 magnitudes of -band extinction). Table 5 lists the results.

One may compare the values obtained here for CFHT/MegaCam with those provided for a G2V star on the stellar isochrone web site of Padova Observatory (http://stev.oapd.inaf.it/cgi-bin/cmd), where the same extinction law is assumed. Their values of A()/A(V) for Megacam are 1.466, 1.167, 0.860, 0.656, 0.500. No comparison value is available through that web interface for WIRCam/Ks.

Note that the foreground extinction towards Virgo is smaller A(V)=0.1. Uncertainties associated with neglecting color terms or using one or the other reference for extinction coefficients will therefore be smaller than 0.5 %.

### Footnotes

1. Based on observations obtained with WIRCam, a joint project of CFHT, Taiwan, Korea, Canada, France, at the Canada-France-Hawaii Telescope (CFHT) which is operated by the National Research Council (NRC) of Canada, the “Institute National des Sciences de l’Univers of the Centre National de la Recherche Scientifique” of France, and the University of Hawaii.
2. slugcomment: To appear in Astrophysical Journal Supplement Series
3. Several designations are used to denote MegaCam filters. For simplicity, we adopt and , since of all filters only is significantly different from SDSS u.
4. Unless stated otherwise, magnitudes quoted in this paper are given in the AB magnitude system. See Appendix A.1 for the numerical conversions between AB and Vega magnitudes for NGVS and NGVS-IR filters.
5. For simplicity, we will indifferently use the letters and to refer to the MegaCam filter.
6. The total exposure time for the corresponding semester in hours is given below the program identifier.
7. Number of observed frames. Note that the given number of frames refers to all science observations. In particular, there are no separate “sky” observations with our pipeline reduction technique (see Section 3.2.3).
8. Number of frames with seeing lower than .
9. Modified Julian date.
10. For a detailed description see the CFHT/WIRCam web pages at http://www.cfht.hawaii.edu.
11. http://www.cfht.hawaii.edu/Instruments/Imaging/WIRCam/
IiwiVersion2Doc.html
12. http://www.astro.queensu.ca/virgo/
13. See Appendix B for a discussion on the small color dependence of the extinction coefficients.
14. http://www.cfht.hawaii.edu/Instruments/WIRCam/ dietWIRCam.html

### References

1. Ames, A. 1930, Annals of Harvard College Observatory, 88, 1
2. Arimoto, N. 1996, From Stars to Galaxies: the Impact of Stellar Physics on Galaxy Evolution, 98, 287
3. Bekki, K., Couch, W. J., Drinkwater, M. J., & Shioya, Y. 2003, MNRAS, 344, 399
4. Bertin, E., Mellier, Y., Radovich, M., et al. 2002, Astronomical Data Analysis Software and Systems XI, 281, 228
5. Bertin, E. 2006, Astronomical Data Analysis Software and Systems XV, 351, 112
6. Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393
7. Bertin, E. 2011, Astronomical Data Analysis Software and Systems XX, 442, 435
8. Binggeli, B., Sandage, A., & Tammann, G. A. 1985, AJ, 90, 1681
9. Bielby, R., Hudelot, P., McCracken, H. J., et al. 2012, A&A, 545, A23
10. Blakeslee, J. P., Jordán, A., Mei, S., et al. 2009, ApJ, 694, 556
11. Blakeslee, J. P. 2012, Ap&SS, 341, 179
12. Cohen, J. G., Blakeslee, J. P., & Ryzhov, A. 1998, ApJ, 496, 808
13. Côté, P., Blakeslee, J. P., Ferrarese, L., et al. 2004, ApJS, 153, 223
14. Cutri, R. M., Skrutskie, M. F., van Dyk, S., et al. 2003, VizieR Online Data Catalog, 2246, 0
15. Dahlen, T., Mobasher, B., Faber, S. M., et al. 2013, ApJ, 775, 93
16. Daddi, E., Cimatti, A., Renzini, A., et al. 2004, ApJ, 617, 746
17. Dalton, G. B., Caldwell, M., Ward, A. K., et al. 2006, Proc. SPIE, 6269
18. Faber, S. M., Willmer, C. N. A., Wolf, C., et al. 2007, ApJ, 665, 265
19. Fellhauer, M., & Kroupa, P. 2002, MNRAS, 330, 642
20. Ferrarese, L., Côté, P., Cuillandre, J.-C., et al. 2012, ApJ, 200, 4
21. Fioc, M., & Rocca-Volmerange, B. 1997, A&A, 326, 950
22. Forte, J. C., Faifer, F. R., Vega, E. I., et al. 2013, MNRAS, 431, 1405
23. Gavazzi, G., Boselli, A., Donati, A., Franzetti, P., & Scodeggio, M. 2003, A&A, 400, 451
24. Hanes, D. A., & Brodie, J. P. 1986, ApJ, 300, 279
25. Hanes, D. A., Côté, P., Bridges, T. J., et al. 2001, ApJ, 559, 812
26. Harris, W. E., Whitmore, B. C., Karakla, D., et al. 2006, ApJ, 636, 90
27. Haşegan, M., Jordán, A., Côté, P., et al. 2005, ApJ, 627, 203
28. Hempel, M., & Kissler-Patig, M. 2004, A&A, 428, 459
29. Hilker, M. 2011, EAS Publications Series, 48, 219
30. Jensen, J. B., Tonry, J. L., Barris, B. J., et al. 2003, ApJ, 583, 712
31. Jones, J. B., Drinkwater, M. J., Jurek, R., et al. 2006, AJ, 131, 312
32. Jordán, A., Côté, P., West, M. J., & Marzke, R. O. 2002, ApJ, 576, L113
33. Jordán, A., Peng, E. W., Blakeslee, J. P., et al. 2009, ApJS, 180, 54
34. Kissler-Patig, M., Brodie, J. P., & Minniti, D. 2002, A&A, 391, 441
35. Kotulla, R., Fritze, U., Weilbacher, P., & Anders, P. 2009, MNRAS, 396, 462
36. Lançon, A., & Mouhcine, M. 2000, Massive Stellar Clusters, 211, 34
37. Lane, K. P., Almaini, O., Foucaud, S., et al. 2007, MNRAS, 379, L25
38. Lawrence, A., Warren, S. J., Almaini, O., et al. 2007, MNRAS, 379, 1599
39. Lawrence, A., Warren, S. J., Almaini, O., et al. 2012, VizieR Online Data Catalog, 2314, 0
40. Lin, L., Dickinson, M., Jian, H.-Y., et al. 2012, ApJ, 756, 71
41. Liu, Y., Zhou, X., Ma, J., et al. 2005, AJ, 129, 2628
42. Ly, C., Malkan, M. A., Kashikawa, N., et al. 2012, ApJ, 757, 63
43. McCracken, H. J., Capak, P., Salvato, M., et al. 2010, ApJ, 708, 202
44. McCracken, H. J., Milvang-Jensen, B., Dunlop, J., et al. 2012, A&A, 544, A156
45. Merson, A. I., Baugh, C. M., Helly, J. C., et al. 2013, MNRAS, 429, 556
46. Mouhcine, M., González, R. A., & Liu, M. C. 2005, MNRAS, 362, 1208
47. Kim, H.-S., Yoon, S.-J., Sohn, S. T., et al. 2013, ApJ, 763, 40
48. King, I. R. 1978, ApJ, 222, 1
49. Kurk, J., Cimatti, A., Daddi, E., et al. 2013, A&A, 549, A63
50. McDonald, M., Courteau, S., Tully, R. B., & Roediger, J. 2011, MNRAS, 414, 2055
51. Mei, S., Blakeslee, J. P., Côté, P., et al. 2007, ApJ, 655, 144
52. Misgeld, I., & Hilker, M. 2011, MNRAS, 414, 3699
53. Park, C., & Choi, Y.-Y. 2005, ApJ, 635, L29
54. Park, H. S., Lee, M. G., & Hwang, H. S. 2012, ApJ, 757, 184
55. Paudel, S., Duc, P.-A., Côté, P., et al. 2013, ApJ, 767, 133
56. Peng, E. W., Jordán, A., Côté, P., et al. 2008, ApJ, 681, 197
57. Peng, E. W., Jordán, A., Blakeslee, J. P., et al. 2009, ApJ, 703, 42
58. Pessev, P. M., Goudfrooij, P., Puzia, T. H., & Chandar, R. 2008, MNRAS, 385, 1535
59. Pota, V., Forbes, D. A.,Romanowsky, A. J., et al. 2013, MNRAS, 428, 389
60. Puget, P., Stadler, E., Doyon, R., et al. 2004, Proc. SPIE, 5492, 978
61. Puzia, T. H., Zepf, S. E., Kissler-Patig, M., et al. 2002, A&A, 391, 453
62. Puzia, T. H., Mobasher, B., & Goudfrooij, P. 2007, AJ, 134, 1337
63. Raimondo, G., Cantiello, M., Brocato, E., & Biscardi, I. 2010, Memorie della Societa Astronomica Italiana Supplementi, 14, 63
64. Rangel, C., Nandra, K., Laird, E. S., & Orange, P. 2013, MNRAS, 428, 3089
65. Reaves, G. 1956, AJ, 61, 69
66. Reaves, G. 1983, ApJS, 53, 375
67. Riffel, R., Ruschel-Dutra, D., Pastoriza, M. G., et al. 2011, MNRAS, 410, 2714
68. Romanowsky, A. J., Strader, J., Brodie, J. P., et al. 2012, ApJ, 748, 29
69. Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103
70. Schuberth, Y., Richtler, T., Hilker, M., et al. 2012, A&A, 544, A115
71. Sick, J., Courteau, S., Cuillandre, J.-C., et al. 2013, arXiv:1303.6290
72. Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163
73. Strader, J., Romanowsky, A. J., Brodie, J. P., et al. 2011, ApJS, 197, 33
74. Tamura, N., Sharples, R. M., Arimoto, N., et al. 2006a, MNRAS, 373, 588
75. Tamura, N., Sharples, R. M., Arimoto, N., et al. 2006b, MNRAS, 373, 601
76. Taylor, E. N., Hopkins, A. M., Baldry, I. K., et al. 2011, MNRAS, 418, 1587
77. Tonry, J. L., Ajhar, E. A., & Luppino, G. A. 1990, AJ, 100, 1416
78. de Vaucouleurs, G., de Vaucouleurs, A., Corwin, H. G., Jr., et al. 1991, Third Reference Catalogue of Bright Galaxies. Volume I: Explanations and references. Volume II: Data for galaxies between 0 and 12. Volume III: Data for galaxies between 12 and 24., by de Vaucouleurs, G.; de Vaucouleurs, A.; Corwin, H. G., Jr.; Buta, R. J.; Paturel, G.; Fouqué, P.. Springer, New York, NY (USA), 1991, 2091 p., ISBN 0-387-97552-7, Price US\$ 198.00. ISBN 3-540-97552-7, Price DM 448.00. ISBN 0-387-97549-7 (Vol. I), ISBN 0-387-97550-0 (Vol. II), ISBN 0-387-97551-9 (Vol. III)
79. van Dokkum, P. G. 2001, PASP, 113, 1420
80. Vogel, H. 1979, Mathematical Biosciences, 44, 179-189
81. Webb, J. J., Harris, W. E., & Sills, A. 2012, ApJ, 759, L39
82. Worthey, G. 1994, ApJS, 95, 107
83. Yoon, S.-J., Sohn, S. T., Lee, S.-Y., et al. 2011, ApJ, 743, 149
84. Yuma, S., Ohta, K., Yabe, K., Kajisawa, M., & Ichikawa, T. 2011, ApJ, 736, 92
85. Yuma, S., Ohta, K., & Yabe, K. 2012, ApJ, 761, 19
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters