THE NEXT GENERATION VIRGO CLUSTER SURVEY (NGVS) XXXI. THE KINEMATICS OF INTRA-CLUSTER GLOBULAR CLUSTERS IN THE CORE OF THE VIRGO CLUSTER

ALESSIA LONGOBARDI,<fnref n="1" /> ERIC W. PENG,<fnref n="2" />,<fnref n="1" /> PATRICK CÔTÉ,<fnref n="3" /> J. CHRISTOPHER MIHOS,<fnref n="4" />

LAURA FERRARESE,<fnref n="3" /> THOMAS H. PUZIA,<fnref n="5" /> ARIANE LANÇON,<fnref n="6" /> HONG-XIN ZHANG,<fnref n="7" />,<fnref n="5" />

ROBERTO P. MUÑOZ,<fnref n="5" /> JOHN P. BLAKESLEE,<fnref n="8" />,<fnref n="9" /> PURAGRA GUHATHAKURTA,<fnref n="10" /> PATRICK R. DURRELL,<fnref n="11" />

RÚBEN SÁNCHEZ-JANSSEN,<fnref n="12" /> ELISA TOLOBA,<fnref n="13" /> ANDRÉS JORDÁN,<fnref n="5" />,<fnref n="14" /> SUSANA EYHERAMENDY,<fnref n="14" />,<fnref n="15" />

JEAN-CHARLES CUILLANDRE,<fnref n="16" /> STEPHEN D. J. GWYN,<fnref n="8" /> ALESSANDRO BOSELLI,<fnref n="17" /> PIERRE-ALAIN DUC,<fnref n="6" />

CHENGZE LIU,<fnref n="18" /> KARLA ALAMO-MARTÍNEZ,<fnref n="5" /> MATHIEU POWALKA,<fnref n="6" /> AND SUNGSOON LIM<fnref n="2" />

(Accepted 07 14, 2018)

ABSTRACT

Intra-cluster (IC) populations are expected to be a natural result of the hierarchical assembly of clusters, yet their low space densities make them difficult to detect and study. We present the first definitive kinematic detection of an IC population of globular clusters (GCs) in the Virgo cluster, around the central galaxy, M87. This study focuses on the Virgo core for which the combination of NGVS photometry and follow-up spectroscopy allows us to reject foreground star contamination and explore GC kinematics over the full Virgo dynamical range. The GC kinematics changes gradually with galactocentric distance, decreasing in mean velocity and increasing in velocity dispersion, eventually becoming indistinguishable from the kinematics of Virgo dwarf galaxies at $R > 320$ kpc. By kinematically tagging M87 halo and intra-cluster GCs we find that 1) the M87 halo has a smaller fraction ($52 \pm 3\%$) of blue clusters with respect to the IC counterpart ($77 \pm 10\%$), 2) the $(g' - r')_0$ vs $(i' - z')_0$ color-color diagrams reveal a galaxy population that

<footnote n="1">¹Kavli Institute for Astronomy and Astrophysics, Peking University, Beijing 100871, China; alongobardi@pku.du.cn; peng@pku.edu.cn</footnote>

<footnote n="2">²Department of Astronomy, Peking University, Beijing 100871, China</footnote>

<footnote n="3">³NRC Herzberg Astronomy and Astrophysics, National Research Council, 5071 West Saanich Road, Victoria, BC V9E 2E7, Canada</footnote>

<footnote n="4">⁴Department of Astronomy, Case Western Reserve University, Cleveland, OH 41106, USA</footnote>

<footnote n="5">⁵Institute of Astrophysics, Pontificia Universidad Católica de Chile, Av. Vicuña Mackenna 4860, 7820436 Macul, Santiago, Chile</footnote>

<footnote n="6">⁶Observatoire astronomique de Strasbourg, Université de Strasbourg, CNRS, UMR 7550, 11 rue de l'Université, F-67000 Strasbourg, France</footnote>

<footnote n="7">⁷CAS Key Laboratory for Research in Galaxies and Cosmology, Department of Astronomy, University of Science and Technology of China, Hefei 230026, China</footnote>

<footnote n="8">⁸NRC Herzberg Astronomy and Astrophysics, 5071 West Saanich Road, Victoria, BC V9E 2E7, Canada</footnote>

<footnote n="9">⁹Gemini Observatory, Casilla 603, La Serena 1700000, Chile</footnote>

<footnote n="10">¹⁰UCO/Lick Observatory, University of California, Santa Cruz, 1156 High Street, Santa Cruz, CA 95064, USA</footnote>

<footnote n="11">¹¹Department of Physics and Astronomy, Youngstown State University, Youngstown, OH 44555, USA</footnote>

<footnote n="12">¹²UK Astronomy Technology Centre, Royal Observatory Edinburgh, Blackford Hill, Edinburgh, EH9 3HJ, UK</footnote>

<footnote n="13">¹³Department of Physics, University of the Pacific, 3601 Pacific Avenue, Stockton, CA 95211, USA</footnote>

<footnote n="14">¹⁴Millennium Institute of Astrophysics, Av. Vicuña Mackenna 4860, 7820436 Macul, Santiago, Chile</footnote>

<footnote n="15">¹⁵Department of Statistics, Faculty of Mathematics Pontificia Universidad Católica de Chile, Av. Vicuña Mackenna 4860, 7820436 Macul, Santiago, Chile</footnote>

<footnote n="16">¹⁶CEA/IRFU/SAp, Laboratoire AIM Paris-Saclay, CNRS/INSU, Université Paris Diderot, Observatoire de Paris, PSL Research University, F-91191 Gif-sur-Yvette Cedex, France</footnote>

<footnote n="17">¹⁷Aix Marseille Université, CNRS, LAM (Laboratoire d'Astrophysique de Marseille) UMR 7326, F-13388 Marseille, France</footnote>

<footnote n="18">¹⁸Department of Astronomy, Shanghai Key Laboratory for Particle Physics and Cosmology, Shanghai Jiao Tong University, Shanghai 200240, China</footnote>
is redder than the IC population that may be due to a different composition in chemical abundance and progenitor mass, and 3) the ICGC distribution is shallower and more extended than the M87 GCs, yet still centrally concentrated. The ICGC specific frequency, $S_{N,ICL} = 10.2 \pm 4.8$, is consistent with what is observed for the population of quenched, low-mass galaxies within 1 Mpc from the cluster's center. The IC population at Virgo's center is thus consistent with being an accreted component from low-mass galaxies tidally stripped or disrupted through interactions, with a total mass of $M_{ICL,tot} = 10.8 \pm 0.1 \times 10^{11} M_{\odot}$.

*Keywords:* galaxies: clusters: individual (Virgo) — galaxies: individual (M87) — galaxies: kinematics — galaxies: star clusters: general — globular clusters: general
# 1. INTRODUCTION

Mergers are believed to play a dominant role in the assembly of massive galaxies (De Lucia & Blaizot 2007). A two-phase formation scenario (Naab et al. 2009) predicts that in a first dissipative phase (high redshift, $z \le 2$) stars form quickly and build the compact innermost regions of massive galaxies (with stars rich in α-elements; Thomas et al. 2005); later, minor mergers dominate the mass assembly in the outermost regions (Oser et al. 2010) and the accretion of metal-poor stars from smaller stellar systems results in the variation of the galaxy stellar properties. In dense environments, accretion is even more dramatic, such that close to the dynamical center of clusters, central galaxies are expected to have accreted a majority of their stars (Cooper et al. 2015).

In galaxy clusters, a fraction of their baryonic content is represented by the intra-cluster light (ICL, or “diffuse light”), a component that is gravitationally unbound to cluster galaxies, but bound to the cluster potential (e.g., Murante et al. 2004; Rudick et al. 2006; Dolag et al. 2010; Puchwein et al. 2010; Contini et al. 2014). Galaxy interactions, as well as tidal interactions between galaxies and the cluster potential, are both believed to play an important role in the production of the ICL. Hence, the ICL is thought to be intimately linked to the dynamical history of the cluster, so that the ICL’s observable features contain information about the evolutionary processes that took place in these dense environments.

Despite its dynamical definition, and low surface brightness (SB) (its peak in SB correspond to ~ 1% of the brightness of the night sky, Vílchez-Gómez 1999), the ICL is usually identified on the basis of its photometric properties. It is either identified as any optical light below a fixed SB limit (Feldmeier et al. 2004a; Mihos et al. 2005; Zibetti et al. 2005; Mihos et al. 2017), or as the less concentrated light profile that overlaps onto the one of the galaxy’s halo (Gonzalez et al. 2005; Seigar et al. 2007). These studies have shown that the detected ICL is often a discernibly separate entity from the host galaxies, with well defined transitions in the surface brightness profile, axis ratio, and position angle, and whose evolution is tied to the cluster as whole rather than to the central galaxy. In theoretical studies, in which the ICL is dynamically identified in terms of either binding energy (stars that are not bound to identified galaxies, including the central one; Murante et al. 2004) or velocity distribution (broader for the ICL than for the galaxy halo; Dolag et al. 2010; Contini et al. 2014), the authors find two distinct stellar populations in terms of kinematics, spatial distribution and physical properties like age and metallicity, suggesting that the identification of the system galaxy+ICL as a single en-

tity is ill-defined. However, in several observational and theoretical analysis (e.g., Gonzalez et al. 2007; Cooper et al. 2015), no separation is made between galaxy halos and ICL, and the two components are treated as a continuum.

In nearby clusters, individual tracers can be used to study the diffuse cluster light around galaxies. They can trace its spatial distribution, age, metallicity, and, when spectroscopy is available, they allow us to gather information on the kinematics of this diffuse component, leading to a less ambiguous definition of galaxy halos and ICL (Longobardi et al. 2015a). Individual stellar tracers include supernovae (Gal-Yam et al. 2003; Dilday et al. 2010; Sand et al. 2011; Barbary et al. 2012), red-giant and super-giant stars (Ferguson et al. 1998; Durrell et al. 2002; Ohyama & Hota 2013), planetary nebulae (PNe; Arnaboldi et al. 1996; Feldmeier et al. 2004a; Arnaboldi et al. 2004; Gerhard et al. 2005; Aguerri et al. 2005; Ventimiglia et al. 2011; Longobardi et al. 2013, 2015a; Hartke et al. 2017, 2018), and globular clusters (GCs; McLaughlin et al. 1994; Jordán et al. 2004; Williams et al. 2007; Lee et al. 2010; Peng et al. 2011; Strader et al. 2011; Romanowsky et al. 2012; Durrell et al. 2014; Alamo-Martínez & Blakeslee 2017; Ko et al. 2017).

GCs are compact groups of stars that are found to inhabit all types of galaxies more luminous than ~ $3 \times 10^6 L_\odot$ (e.g., Georgiev et al. 2010), with the most luminous systems characterized by higher GC specific frequencies (e.g., Harris & van den Bergh 1981; Brodie & Strader 2006; Peng et al. 2008; Georgiev et al. 2010). Several GC studies have shown that the properties of these systems are correlated with those of their host galaxies (Harris 1991; Brodie & Huchra 1991; Peng et al. 2006, 2008; Harris et al. 2013), suggesting a link between their formation and the evolution of the galaxy itself. It is now well established that the GC systems, within the majority of large galaxies, possess two or more subpopulations of clusters characterized by very different chemical composition (e.g., Gebhardt & Kissler-Patig 1999; Puzia et al. 1999; Forbes & Forte 2001; Kundu & Whitmore 2001; Larsen et al. 2001; Puzia et al. 2005a,b; Peng et al. 2006; Leaman et al. 2013; Tonini 2013; Harris et al. 2017). The most popular division uses optical colors to identify a red (metal-rich) and a blue (metal-poor) population of GCs, although the existence of an additional, intermediate, population of GCs has been claimed (Peng et al. 2009; Strader et al. 2011; Agnello et al. 2014). To these different classes correspond different observed properties, so that red and blue GCs are measured to have different kinematics and spatial distribution (Zepf et al. 2000; Côté et al. 2001; Schroder
et al. 2002; Côté et al. 2003; Perrett et al. 2003; Peng et al. 2004; Schuberth et al. 2010; Strader et al. 2011; Coccato et al. 2013; Agnello et al. 2014; Zhang et al. 2015). Also, the relative fraction of blue-to-red GCs increases with galactocentric radius, probably tracing the hierarchical processes that built up the galaxies stellar halos (e.g., Côté et al. 1998, 2000, 2002; Tonini 2013; Lim et al. 2017). It is important to emphasize, however, that the interpretation of a single color as a measure of metallicity is a simplification. Such an interpretation presents the uncertainty driven by the nonlinearity in the relation (Peng et al. 2006; Cantiello & Blakeslee 2007; Richtler 2013), and wrongly assumes that the same stellar population characteristics found in the Local Group apply to dense environments (Blakeslee et al. 2012; Powalka et al. 2017).

As consequence of their brightness and spatial extension, GCs are powerful tracers of the galaxies' outer halos, also bridging the transition region between the galaxies and the ICL that surrounds them. Single intra-group/cluster GCs have been detected in A1185 (Jordán et al. 2003; West et al. 2011), Fornax (Firth et al. 2008; Schuberth et al. 2008), Virgo (Williams et al. 2007; Firth et al. 2008; Lee et al. 2010; Ko et al. 2017), Coma (Peng et al. 2011), and A1689 (Alamo-Martínez & Blakeslee 2017). However, the identification of a large sample of tracers has been elusive, mainly due to a significant level of contamination (see Sect. 2.3).

The Virgo cluster, the nearest large-scale structure in the Universe, and its central galaxy M87 have long been the targets of GC studies with the aim of tracing their formation and evolution (e.g., Baum 1955; Mould et al. 1987; Jordán et al. 2002; Côté et al. 2004; Ferrarese et al. 2012; Brodie et al. 2014). Virgo is characterized by both spatial and kinematic substructures (Binggeli et al. 1987). Moreover, the evidence that many galaxies are presently falling towards the cluster core (Tully & Shaya 1984; Conselice et al. 2001; Boselli et al. 2008, 2014), as well as the presence of a complex network of extended tidal features (Mihos et al. 2005; Mihos et al. 2017), suggest that Virgo core is not completely in equilibrium. Close to its dynamical center lies M87 (Binggeli et al. 1987; Nulsen & Bohringer 1995; Mei et al. 2007). It is considered a cD galaxy (e.g., Carter & Dixon 1978) with an extended stellar envelope that reaches out to projected galactocentric radii of ~150 kpc (Ferrarese et al. 2006; Kormendy et al. 2009; Janowiecki et al. 2010). Evidence of the galaxy experiencing a gas-rich major-merger event can be found in its central regions, hosting a kinematically distinct core (KDC) in its inner 5 kpc (Emsellem et al. 2014). Moreover Hα studies have shown the presence of prominent filaments of ionized gas

extending out to 10 kpc from the galaxy's nucleus (Arp 1967; Ford & Butcher 1979; Sparks et al. 1993; Gavazzi et al. 2000). Evidence of accretion are traced by age and metallicity gradients in the galaxy's inner region (Liu et al. 2005; Montes et al. 2014), together with blue color gradients towards the outer regions (Rudick et al. 2010; Mihos et al. 2017). Furthermore, at large radii, kinematic signatures of accretion events are found in the orbital distribution of GCs (Agnello et al. 2014; Oldham & Evans 2016), in the kinematics of UCDs (Zhang et al. 2015), as well as in the presence of kinematic substructures in the velocity phase-space of GCs (Romanowsky et al. 2012) and PNs (Longobardi et al. 2015b). A large accretion event (the crown of M87, Weil et al. 1997; Longobardi et al. 2015b) was found to contribute 60% of the total light at its highest density point (60–90 kpc NW of the galaxy; Longobardi et al. 2015b), showing that M87 is still assembling in a substantial way.

This work follows other recent Next Generation Virgo Cluster Survey analyses with the aim of studying the Virgo core (Zhu et al. 2014; Grossauer et al. 2015; Liu et al. 2015a; Zhang et al. 2015; Sánchez-Janssen et al. 2016; Ferrarese et al. 2016; Roediger et al. 2017), and the correlations between the GC properties and the evolution of the environment that surround them (e.g., Powalka et al. 2016a,b, 2017).

We make use of the deep and extended photometric information provided by the Next Generation Virgo Cluster surveys (NGVS/NGVS-IR Ferrarese et al. 2012; Muñoz et al. 2014) in the central 2x2 deg² of Virgo, to study the transition region between the central galaxy and the IC space. The aim is to identify and separate the GC populations associated with these components and study their properties separately. The paper is structured as follows: in Section 2 we present our data, the different sources of contamination, hence our working sample of GCs. In Section 3 we study the GC line-of-sight velocity distribution (LOSVD) and kinematically separate the M87 halo component from the Virgo ICL. The properties of both M87 halo and Virgo intra-cluster GCs are analyzed in Section 4 in terms of colors and spatial distributions. We discuss our results in Section 5 and finally report our conclusions in Section 6.

Throughout the paper we assume a distance for Virgo and M87 of 16.5 Mpc (Mei et al. 2007; Blakeslee et al. 2009), implying a physical scale of 80 pc arcsec⁻¹.

## 2. OBSERVATIONS & DATA

This work examines the GC population in the core region of the Virgo cluster through a synergy between optical/near-infrared photometry from the Next Generation Virgo Cluster Survey (NGVS, Ferrarese et al.
2012)/ Next Generation Virgo Cluster Survey-Infrared (NGVS-IR, Muñoz et al. 2014), and spectroscopy, the latter coming from a compilation of different campaigns. In what follows we give a brief description of these surveys and refer the reader to the references therein for additional details.

## 2.1. NGVS photometry

The NGVS is an optical imaging survey of the Virgo cluster carried out with MegaCam on the CFHT (Ferrarese et al. 2012). It covers a total area of 104 deg$^2$, covering the two main sub-clusters (Virgo A centered on M87, and Virgo B centered on M49) out to their virial radii. Ferrarese et al. (2012) adopted several methods of background subtraction and image combination to produce stacked NGVS images. Here we work with the stacks built using the MegaPipe global background subtraction, combined with the artificial skepticism algorithm (Stetson et al. 1989). This provides high accuracy photometry for sources of small spatial extent, making the NGVS a deep photometric survey in the $u^*g'r'i'z'$ bands, that for point sources reaches a depth of $g' = 25.7$ mag$^1$. In addition, images in all filters have a FWHM that is lower than 1.0, with the best seeing conditions obtained for the $i'$ band with an upper limit of 0.6. All final images have the same astrometric reference frame, tied to the positions of stars in the Sloan Digital Sky Survey, and the same spatial resolution with a physical scale of 186 px$^{-1}$.

For the Virgo core, a 3.62 deg$^2$ region centered on M87, a deep $K_s$ band imaging survey (NGVS-IR) was also carried out with WIRCam on the CFHT (Muñoz et al. 2014). Any raw images with a seeing worse than 0.7 were rejected before stacking. This made it possible to produce stacked images with the same pixel scale as the MegaCam stacks, although the original WIRCam pixel scale is 0.3 px$^{-1}$. The Swarp software (Bertin et al. 2002) coupled with a Lanczos-2 interpolation method was used to stack the sky-subtracted images. Thus, the NGVS-IR has a limiting magnitude of $K_s \sim 24.4$ AB mag for point sources, and an image quality better than 0.7. Photometry of all objects on the final, stacked NGVS/NGVS-IR images was conducted using Source Extractor (Bertin & Arnouts 1996). The aperture corrections (due to small variations in the stellar point-spread function across the field) were derived using samples of bright, unsaturated stars in each field

for each filter (more details can be found in Durrell et al. 2014; Liu et al. 2015a; Muñoz et al. 2014).

## 2.2. Spectroscopy

We have compiled a catalog containing all the known redshifts within the NGVS footprint. It combines data from the literature (e.g., Binggeli et al. 1985; Hanes et al. 2001; Strader et al. 2011), public archives such as SDSS (DR12, Alam et al. 2015) and NED$^2$, as well as data from our own spectroscopic surveys dedicated to the detection of compact objects in Virgo. In this work we only focus on the inner $2^\circ \times 2^\circ$ around M87, and in what follows we will give a brief description of the surveys that provided data for this analysis.

### Spectroscopic data in the Virgo core—

**The Literature sample.** In addition to our own surveys (described below), the central region of the Virgo cluster has been the target of several spectroscopic campaigns. Strader et al. (2011) (S11) combined the pre-2011 published data together with their new observations to compile a sample of radial velocities covering M87 within a 40' radius. As a compilation of different surveys (e.g., Hanes et al. 2001; Côté et al. 2001), data from the S11 catalog come from observations carried out with different instruments, hence characterized by different resolutions and spectral coverage. The authors extensively describe the different observational strategies and the compilation procedure that led to the creation of the S11 catalog, hence we refer the reader to their work for a detailed description of the surveys. Here we note that we use the radial velocities provided in their catalog, but in some cases have modified the object classification based on our high-quality photometric information (see Sec. 2.3.1).

**The MMT campaign.** In 2009 and 2010 we used the Hectospec multifiber spectrograph on the 6.5m MMT telescope to survey the compact stellar systems in the central $2^\circ \times 2^\circ$ area around M87 (Zhang et al. 2015). Data covered the spectral region within 3650Å–9200Å with a resolution of $R = 1000$. By exposing for ~2 hours, the survey depth is $g' \sim 22.5$ mag and velocities have a typical uncertainty of $\sigma_{V_{LOS}} = 30$ km s$^{-1}$.

<footnote>$^1$ Deep imaging in the $u^*$ and $r'$ bands do not cover the entire NGVS footprint, with the $r'$ band that is limited to the central 2 by 2 square degree around M87. Magnitudes are on the CFHT MegaPrime photometric system.</footnote>

<footnote>$^2$ The NASA/IPAC Extragalactic Database (NED) is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.</footnote>
*The AAT campaign.* Compact objects in Virgo were further followed up spectroscopically in March/April 2012 when we surveyed the central region of the Virgo cluster (Virgo A sub-cluster), using the 2dF AAOmega multi-fiber spectrograph on the AAT. Spectra were obtained in the wavelength range 3700Å–8000Å and with a spectral resolution of R = 1300 (Zhang et al. 2015). The survey consisted of nine 2dF pointings with a typical exposure time of 1.5 hours, covering a total sky area of ~30 deg$^2$. At the limiting magnitude of $g' \sim 20.5$ mag and with such an observational set-up the typical uncertainty on the velocity estimates is $\sigma_{V_{LOS}} = 30 \text{ km s}^{-1}$.

*The Keck campaign.* In April 2013 additional data were collected using the DEIMOS spectrograph located at the Keck II 10 m telescope. The instrumental configuration provided a wavelength coverage of 4800Å–9500Å with a spectral resolution of 2.8Å (FWHM). This campaign provided velocity measurements with a mean precision of $\sigma_{V_{LOS}} = 10 \text{ km s}^{-1}$ and down to $g' \sim 24.7$ mag (Toloba et al. 2016). 24.7 must be the faintest object, but certainly not representative of the GC population.

*The Magellan campaign.* Finally, in March 2016 we used the IMACS multi-slit spectrograph on the Magellan Baade telescope to obtain spectra for compact objects in Virgo with $g' < 22$ mag. With a field of view of $27.5' \times 27.5'$, the observations surveyed two regions, one centered on M87 and the second one offset by 30' towards the NW along the major axis of the galaxy. Data were obtained in the wavelength range 3900Å – 9000Å with a spectral resolution of 6.5Å. With a total integration time of 3.5 hrs this survey provided data with a mean velocity uncertainty $\sigma_{V_{LOS}} = 30 \text{ km s}^{-1}$ (Zhang et al. 2018).

As we show below the NGVS photometry provides vital information needed to select the GC sample that we use to analyze the kinematics of the GC system and identify the LOSVDs of the different dynamical components. Hence our working sample is selected from the matched objects between the NGVS photometric and spectroscopic catalogs in the pilot region, resulting in a total sample of $N_{obj} = 2809$ sources. We further decide to analyze objects down to $g' \le 24.0$ mag and $\sigma_{V_{LOS}} \le 50 \text{ km s}^{-1}$. This leads to a total number of sources $N_{obj} = 2551$, with median photometric uncertainties $\sigma_{mag} = [0.01, 0.005, 0.004, 0.005, 0.008, 0.02]$ in the $u^*g'r'i'z'$ and $K_s$ bands, respectively. Of these

objects, 897 have high photometric probability to be a GC (see Sec. 2.3.1). Among these, 592 were gathered by S11, while 232, 9, 29, and 18 radial velocities were observed for the first time in the MMT, AAT, Keck, and Magellan campaign, respectively (the remaining velocities are from public archives). We emphasize that from repeated measurements in different campaigns the velocity estimates are in agreement within the spectroscopic uncertainty threshold of $\sim \sigma_{V_{LOS}} = 50 \text{ km s}^{-1}$. The two spectroscopic catalogs with the highest number of GCs in common (67 common objects from the S11 and MMT samples) only present two velocity estimates that deviate $\sim \sim 1.5 \sigma_{V_{LOS}}$.

## 2.3. Objects classification

The selection of a bona-fide sample of GCs is a crucial first step for our kinematic analysis. The two main types of contaminants in our sample are i) background galaxies and ii) foreground Milky Way stars. A third, more subtle, type of potential contaminants are iii) Virgo UCD galaxies. In what follows, we examine such a contamination and its contribution to the sample of GCs.

### 2.3.1. Background galaxies/Foreground stars contamination

Muñoz et al. (2014) showed that background galaxies, GCs, and foreground stars define different regions in the $u^*i'K_s$ color-color diagram. In Fig. 1, this relation is shown for the photometric sample in the NGVS pilot region. From redder to bluer ($i' - K_s$) colors, we can see the different contribution from background galaxies with various star forming histories at redshifts up to $z \simeq 1$, GCs (which merge into the redshift sequence of passive galaxies at the red end), and foreground main sequence stars. However, a successful compilation of a GC sample implies accounting for the shared contribution between different populations in the transition regions of these three sequences. To do so, we have used the full photometric information gathered by the NGVS/NGVS-IR surveys, and adapted the "extreme deconvolution" (XD) algorithm from Bovy et al. (2011) to classify our data according to multi-dimensional color and concentration information. The details of this classification procedure will be presented in a future work. Here we give a brief description of the procedure and validate the results in the next Section, with the spectroscopic information we have.

In the pilot region the XD approach uses multi-dimensional Gaussians to model the density of objects in the 5-dimensional space: $(u^*-g')_0$, $(g'-i')_0$, $(i'-z')_0$, $(i'-K_s)_0$, and a measure of concentration $iC = i'_4 - i'_8$, where $i'_4$ and $i'_8$ are the point source aperture-corrected
<figure><img src="image_1.png" /><figcaption>**Figure 1.** *u*\*i\'Ks color-color diagram for the NGVS photometric sample in the pilot region with *g*' ≤ 24.0 AB mag. In this plane background galaxies, GCs, and foreground MW stars separate into different regions (Muñoz et al. 2014). Magenta dots identify objects with high photometric probability, *p*<sub>gc</sub> > 0.5, to be GCs (see text for more details).</figcaption></figure>

i-band magnitudes measured in fixed apertures of di-
ameter *r* = 4 pixels (0.''72) and *r* = 8 pixels (1.''44),
respectively. The i-band seeing in the NGVS images
has FWHM < 0.''6 with a median seeing of 0.''54. For
point-like objects *iC* = 0, and extended objects have
*iC* > 0. The Gaussians in the model are convolved with
the observational uncertainties of the data. Therefore,
each data point is assumed to be drawn from the model
convolved with its own set of measurement uncertain-
ties, and at the end of the XD procedure we obtain for
each object the probabilities, *p*<sub>galx</sub>, *p*<sub>gc</sub>, and *p*<sub>star</sub>, to be
either a galaxy, a GC or a star, respectively. Magenta
dots in Fig. 1 are objects with high photometric proba-
bility, *p*<sub>gc</sub> > 0.5, to be GCs.

2.3.2. Spectroscopic validation of the photometric selection

**Background galaxies**—There is a well-defined gap be-
tween the Virgo members and background galaxies at ra-
dial velocities $V_{\text{galx}} = 3000 \text{ km s}^{-1}$ (Binggeli et al. 1987;
Conselice et al. 2001; Kim et al. 2014). Therefore, the
background component is identified by selecting objects
with $V_{\text{LOS}} > 3000 \text{ km s}^{-1}$. We find that 721 sources in
our sample satisfy this criterion, representing ~ 28%
of our working sample. The XD photometric classifica-
tion identifies these contaminants successfully: > 99%

of the background objects are assigned high photometric
probability to be background galaxies, with mean value
⟨p<sub>galx</sub>⟩ = 0.98, and only 1% would result in having high
photometric probability to be GCs.

**Foreground stars**—The velocity distribution of the Galactic stellar halo is measured to peak at a heliocentric radial velocity of ~40 km s⁻¹ with a velocity dispersion of ~120 km s⁻¹ (e.g., Battaglia et al. 2005). Therefore, the MW star LOSVD overlaps with the LOSVD of objects in Virgo, the latter measured to cover a velocity range –500 km s⁻¹ ≤ V<sub>LOS</sub> ≤ 3000 km s⁻¹ (Binggeli et al. 1987; Conselice et al. 2001). The contamination by foreground stars is visible as a strong and narrow Gaussian that peaks at around 40 km s⁻¹ in Fig. 2 (left-hand panel), where we plot the LOSVD of our working sample with V<sub>LOS</sub> ≤ 3000 km s⁻¹.

To isolate this kinematic component, we run a Gaussian Mixture Model (GMM), i.e., a probabilistic model that divides the observed LOSVD into a linear combination of *K* independent Gaussian probability density functions. At the end of the procedure, the posterior probability, Γ, for a data value to belong to either of the *K* Gaussian components are returned (Fig. 2 right-hand panel, and see Sec.3 for details on the GMM). Tests with the Bayesian Information Criterion (BIC) and Akaike information criterion (AIC) show that the GMM best-fit model is a mixture of three Gaussians (Fig. 2 central panel), identifying the Galactic stellar contamination as component with $V_{\text{star}} = 41 \pm 4 \text{ km s}^{-1}$ and $\sigma_{\text{star}} = 114 \pm 3 \text{ km s}^{-1}$. In Sec. 3 we will show that the remaining two Gaussians are tracing the Virgo cluster and the M87 halo distributions, however as for the moment our goal is to validate the photometric classification, we only isolate the objects with velocities that fall within the stars’ LOSVD, and consider as one single class the remaining objects.

In Fig. 3 the (*u*\* − *i*') vs (*i*′ − *K*<sub>s</sub>) diagram encodes
the spectroscopic information in the color: objects with
higher GMM posterior probabilities to belong to the star
LOSVD are shown in yellow and orange. From this plot
we can see that the majority of these objects are con-
sistent with being foreground stars on the basis of their
*u*\*i'K<sub>s</sub> properties (yellow circles). However, a fraction
of objects with *u*\*i'K<sub>s</sub> colors consistent with being GCs
have velocities that fall within the stars LOSVD. Their
total number is 32, i.e., ~4% of the sample of star can-
didates based on their kinematics. In the bottom panel
of Fig. 3, where the color-color diagram is re-plotted
by weighting the contribution of each data point by the
XD probability, p<sub>gc</sub>, to be a GC, we show that we can
successfully retrieve this fraction using the photomet-
ric classification previously described. Going back to
<figure><img src="image_2.png" /></figure>

**Figure 2.** Left Panel: LOSVD of the spectroscopically confirmed objects in the NGVS pilot region that define our sample. The GMM best-fit model describes the total LOSVD as a mixture of three Gaussians (black continues line) centered on the star, Virgo ICL, and M87 halo distributions, respectively (black dashed lines). Central Panel: BIC and AIC scores for the Gaussian mixtures as function of the number of components. Both the information criteria prefer three mixtures to describe the data as in the left panel. Right Panel: Posterior probabilities that an object is drawn from one of the three components, i.e. the stars (light gray, class 1), the Virgo IC (gray, class 2), and the M87 galaxy (dark gray, class 3) distributions, as function of its LOSV. For a given velocity value, the vertical extent of each region is proportional to that probability.

the top panel of Fig. 3, we also see that there are few misclassified objects (~9% of our working sample with $V_{LOS} \le 3000 \text{ km s}^{-1}$) with high photometric probability, $p_{galx}$, to be background galaxies (cyan and orange triangles). These objects have been visually inspected and reintroduced in the bona-fide sample of GCs when their velocity is consistent with the Virgo IC/M87 halo components (cyan triangles, 85 objects). For those whose velocities fall within the star LOSVD (orange triangles, 77 objects), we only reintroduce those that are in the locus of the color-color diagram where the majority of GCs sits (4 objects). All the reintroduced objects have been assigned $p_{gc} = 1.0$. To summarize, we have shown that the photometric information gathered by the NGVS is enough to properly separate foreground stars and Virgo compact objects, allowing us to avoid hard cuts in the velocity distribution for $V_{LOS} < 500 \text{ km s}^{-1}$.

### 2.3.3. Virgo UCD contamination

In the last 20 years evidence has emerged for the existence of a family of compact stellar systems, called “ultra-compact dwarfs” (UCDs). Their sizes and luminosities are intermediate between compact elliptical galaxies and the most massive GCs (Hilker et al. 1999; Drinkwater et al. 2000; Haşegan et al. 2005; Misgeld &

Hilker 2011), and their colors and stellar populations are similar to GCs and nuclear star clusters (Taylor et al. 2010; Spengler et al. 2017). Recently, Zhang et al. (2015) studied the properties of a sample of UCD galaxies within ~1 deg of M87, and compared it to the properties of the red and blue GC populations. They found that they do not share the same kinematics as the GC system, however no distinction was made between galaxy halo and intra-cluster GCs, leading to the question as to whether UCDs may trace the IC component instead, hence sharing similar properties as the ICGCs. The evidence that the Virgo UCDs never reach the cluster dispersion, nor do they trace the Virgo systemic velocity, is enough to argue that these systems are not to be associated with the ICL. However, the argument gets stronger if we consider that at large radii the UCDs behave more closely to the red GCs (in terms of density profiles and rotational axis). As we will show in this work (see Sec. 3), the presence of an IC component particularly affects the GC properties at large radii and it is dominated by the blue GC population (see Sec. 4). Thus, if the UCDs did trace the ICL they would rather show similar density distributions and kinematics as the blue GC population at large radii, that we do not observe.
<figure><img src="image_3.png" /></figure>

**Figure 3.** Top Panel: $u^*i'K_s$ diagram color coded based on the spectroscopic information: yellow dots are objects with higher posterior probabilities to belong to the star LOSVD (class 1), while black crosses are more likely to be drawn from the M87 galaxy/ICL distributions (class 2/3). Magenta triangles represent background galaxies. Cyan and orange triangles are objects with velocities consistent with being part of class 1 (orange) or class 2/3 (cyan), however they have photometric probability, to be background galaxies. Bottom Panel: Same as Top panel with the objects size scaled by its photometric probability to be a GC. The photometric information allows us to retrieve the fraction of objects with colors consistent with being GCs but with velocities within the stars LOSVD (yellow dots). Only 1% of background galaxies (magenta triangles) is misclassified as GCs following the photometric classification.

We then use the information gathered by Zhang et al. (2015); Liu et al. (2015a) to exclude the UCDs from our sample of bona-fide GCs. These authors identified UCDs based on their photometric properties, such as size and colors. By cross matching the catalogs from Zhang et al. (2015) and Liu et al. (2015a) with our sample of compact objects we identify 85 UCD candidates that we remove from our sample of bona-fide GCs. The selection of these objects is reliable, only missing ~6% of genuine UCDs for $g' < 21.5$ (Zhang et al. 2015; Liu et al. 2015a). However, under the assumption that UCDs are the surviving nuclei of tidally stripped galaxies, Ferrarese et al. (2016) argued that we expect a contribution

in the $21.5 < g' \le 24.0$ mag range that amounts to ~100 objects. Given our spectroscopic incompleteness (on average we only observe 10% of the photometric candidates, see Sec. 4.2), this implies a residual contamination by UCDs of ~10 objects in our spectroscopic sample. Such a fraction is negligible with respect to the final sample of GCs, and does not affect any of the scientific results that we present in this study.

### 2.3.4. GCs bound to other Virgo members

In the previous sections we have subtracted the contribution of stars, background galaxies and Virgo UCDs from our working sample. However, as our goal is to study the GC population of the M87 halo and the Virgo IC component we also have to flag GCs that are likely to be bound to other Virgo members. To do this we considered GCs to belong to another galaxy if 1) The GC velocity relative to the galaxy velocity, $V_g$, is within $3 \times \sigma_g$, where $\sigma_g$ is the velocity dispersion of the galaxy (assuming a Gaussian LOSVD, 99.7% of the distribution falls within this limit), and 2) their distance from the galaxy's center is within 10 effective radii, $R_e$, of the galaxy. This generous radial criterion was chosen to account for possible differences between the galaxy light and the GC spatial distribution. Fits to the surface brightness and surface density profiles of galaxies' light and GC systems show that the latter are, on average, more extended than that of the host galaxies' (e.g., Puzia et al. 2004; Kartha et al. 2014; Forbes 2017; Hudson & Robison 2018). For less massive, dwarf galaxies the scaling relation can be reduced to $R_{e,GC} = 1.5 \times R_{e,gal}$ (e.g., Peng & Lim 2016).

Within a 2° radius around M87 there are 197 Virgo galaxies with known velocity, 129 of which are early-type dwarf galaxies, as compiled by the GOLDMine project (but see also Toloba et al. (2014), Gavazzi et al. 2003). The GOLDMine catalog also provided us the velocity dispersion for a fraction of galaxies, while the effective radius information was retrieved from McDonald et al. (2011). For those galaxies with no measured velocity dispersion and/or effective radius, mostly dwarfs, we have given a fixed value of $\sigma_g = 50 \text{ km s}^{-1}$ and $R_e = 30''$. At the end of this procedure we have identified 56 GCs as belonging to other Virgo galaxies. Finally, we are aware that there will be a fraction of galaxies in the surveyed area with no velocity information, hence that is not considered in this analysis. However, we note that additional GCs bound to other galaxies would show a correlation in position and velocity that we checked for. This test showed no evidence of any possible spatial/kinematic correlation, thus we can state that if there is a contribution from GCs belonging to other nearby
<figure><img src="image_4.png" /></figure>

**Figure 4.** Spatial distribution of GCs in Virgo with spectroscopic information (purple dots) superimposed on to the g-band NGVS image of the Virgo core. Blue triangles show the position on the sky of UCDs, while orange crosses show GCs that we associated to other Virgo galaxies. Aside M87 (salmon star), we highlight the position of other four main galaxies as given in the plot. Red crosses identify Virgo galaxies with negative $V_{LOS}$. The gray ellipses, oriented as the M87's isophotes (major and minor-axes inclination as depicted by the dotted lines), divide the sky in four regions. North is up, East is to the left.

galaxies in our bona-fide sample of objects, it is negligible.

To summarize, in this section we have analyzed and removed the contribution from contaminants in our working sample to define our final sample of bona-fide GCs. We have found that:

1. Background galaxies (i.e. no Virgo members with $V_{LOS} > 3000 \text{ km s}^{-1}$) represent the ~28% of our working sample;

2. Once the background contamination is subtracted, foreground MW stars contribute ~46% of the data, leaving 975 sources with high photometric probability, $p_{gc}$ to be a compact Virgo object;

3. Among the compact objects, 85, i.e. 9%, are identified as UCDs;

4. Finally, 56 GCs (6% of the total GC sample) are consistent with belonging to other Virgo members.

This selection led to a final sample of 837 confirmed GCs that we will analyze in the next Sections. Their spatial distribution is shown in Fig. 4 (purple dots) together with the position on the sky of UCDs (blue triangles) and GCs associated to other Virgo members (orange crosses; the centers of other Virgo galaxies are signaled by red crosses).

### 3. THE M87 GALAXY HALO AND VIRGO ICL AS TRACED BY GLOBULAR CLUSTERS

When studying the light at the center of a cluster reaching out ~300 kpc (1°) in radius from the central galaxy we are tracing the transition region between the central galaxy halo and the IC component. Longobardi et al. (2015a) analyzed a 0.5 deg² area around M87 using PN kinematics to dynamically identify and separate the galaxy halo and the ICL. They showed that the two components overlap at all radii, but with different PN properties, hence consistent with having different parent stellar populations.

In what follows, we will investigate whether GCs trace a similar galaxy halo-IC interface. If there are, indeed, ICGCs, in what systems did they form and do they show photometric properties that differ from the M87 halo GCs?

#### 3.1. Projected GC velocity phase-space

In Fig. 5 we show the velocity-position phase space diagram (central and bottom panels) for our GC sample. The GC LOS velocities are plotted against the major axis distance, the latter computed as $R = [(x/(1-e))^2 + y^2]^{1/2}$ for a system aligned along the y-axis, with position angle P.A. = 155.0°, and ellipticity $e = 0.4$ (Ferrarese et al. 2006; Janowiecki et al. 2010). The relations for positive and negative R are plotted separately to trace the North-West (NW) (central panel) and South-East (SE) (bottom panel) halves with respect to the center of M87, respectively (see Fig. 4 for a representation of the spatial configuration).

This plot shows a system dominated by random motion, with the GC kinematics concentrated around the M87 systemic velocity $V_{sys} = 1308 \text{ km s}^{-1}$, together with a scattered distribution of high- and low-velocity objects on both sides of the galaxy. As shown in the histograms of the velocities (side panels), the wings (or tails) of the GC LOSVDs are strong and quite asymmetric with respect to a Gaussian distribution (red curve) centered on the M87 systemic velocity and with
<figure><img src="image_5.png" /></figure>

**Figure 5.** **Top Panel:** LOS velocity dispersion profile for the GCs in the NGVS pilot region (dots with the green shaded area delimiting the 1 σ region) and for the dwarf early-type galaxies within a 2° radius from the cluster's center (purple crosses with blue shaded area delimiting the 1 σ region). The GC velocity dispersion profile is observed to increase on both the NW (dark dots and dark shaded area) and SE (yellow dots and light shaded area) sides of M87 reaching the velocity dispersion of the dwarfs in Virgo. **Central Panel:** Projected velocity phase-space that compares the velocity distributions of GCs (black dots) and Virgo dwarf early-type galaxies (purple crosses) in the NW region. The continuous lines trace the running average computed from the total GC sample (green) and the dwarf population (blue). The histogram of the velocities, plotted on the side, shows a LOSVD with strong and asymmetric tails compared to a Gaussian distribution centered on the M87 systemic velocity with a velocity dispersion of ~300 kms⁻¹, and normalized to the total number of GCs (red curve). **Bottom Panel:** Same as central panel, however the relation is shown for the SE region with respect to M87. On both sides the GC mean velocity decreases as function of the distances, deviating from the M87 systemic velocity (Vsys = 1308 kms⁻¹, red line) and merging into the mean velocity of the Virgo dwarfs at |R| ~ 320 kpc. The vertical orange lines separate the phase-space in four elliptical bins, as given in the text.

a velocity dispersion of $\sigma \sim 300 \text{ kms}^{-1}$<fnref n="3" />. Such behavior contrasts with what is observed for non-rotating early-type galaxies, whose LOSVD are symmetric (i.e. $h_3 \sim 0$), hence well described by a single Gaussian to

within few percent (Gerhard 1993; Bender et al. 1994). Orbital anisotropy would create a symmetric variation around the mean velocity (i.e. $h_4 \neq 0$), and is unlikely to be measured at all radii (see Fig. 6 for a plot of the GC LOSVDs as function of the distance from the M87's center). On the other hand, strong tails are seen in simulations, observed to characterize the LOSVD of galaxies in dense environments, and consistent with the

<footnote n="3">³ Here we consider as velocity dispersion representative of M87 the one traced by the red GC population over a range 0.5-30.0'(Zhang et al. 2015)</footnote>
presence of an intra-cluster/group component overlapping with the galaxy halo (Dolag et al. 2010; Longobardi et al. 2015a; Barbosa et al. 2018; Veale et al. 2017). In fact, the relative importance of these “tails” increases as function of the major axis distance, with the GCs gradually deviating from the central galaxy’s kinematics and merging into the cluster potential. This is shown in the same plot by comparing the GC kinematics with the kinematics of the population of early-type dwarfs in Virgo with measured radial velocities (129 objects, Gavazzi et al. 2003). We study the running average of the LOS velocities of GCs (continuous green lines in the central and bottom panels) and early-type dwarfs (continuous blue lines), i.e., we plot the running mean in bins of n objects, with n = 30 for the GCs, n = 35 for the NW dwarfs, and n = 25 for the SE dwarfs. By comparing the different running average profiles it is clear that the GC mean velocity decreases going towards larger distances, deviating more than 2.5 σ at large distances and on both sides from the M87 systemic velocity (red line) and merging into the mean velocity of the Virgo dwarfs at |R| ~ 320 kpc. At these distances the GC velocity dispersion reaches its maximum, σ<sub>LOS</sub> = 590 ± 80 km s<sup>-1</sup>, merging with the cluster potential (Fig. 5, top panel). The comparison between the NW and SE fields around M87 shows that the NW has a more dE galaxies with negative velocities (red crosses in Fig. 4). It is likely that a substantial fraction of these objects comes from the accretion of material that is falling towards the cluster’s center dragged along by the well known infalling M86 subgroup, and it is contributing to what we see today as an IC component.

## 3.2. Kinematic separation between M87 halo and ICL

The comparison with the Virgo dE/dS0 galaxies showed that the GC system around M87 is tracing the transition toward the Virgo cluster potential. Hence, we now further analyze the GC LOSVD, and we do it i) separately for the NW and SE halves, and ii) in four different elliptical bins, or stripes in phase-space (gray ellipses in Fig. 4, or orange lines in Fig. 5), covering the entire GC range. The four regions in space are defined as |R| ≤ 32 kpc, 32 < |R| ≤ 80 kpc, 80 < |R| ≤ 216 kpc, and |R| > 216 kpc, with bin size increasing towards the outer regions to ensure a statistically significant sample of tracers even at large radii, where the GC number density is the lowest. Thus, our choice of radial binning represents the best compromise between similar radial extension and statistical noise. We emphasize that the stability of the results was tested against different choices of radial limits: within the uncertainties the results are consistent with each other. In our analysis,

we assume that both the M87 halo and Virgo IC LOSVDs are Gaussians. As argued in the previous Section, this is a reliable approximation for the galaxy’s stellar velocity distribution, and, from hydrodynamical simulations (Dolag et al. 2010; Cui et al. 2014), Gaussian shapes also describe well the IC component. We are aware, however, that a Gaussian may be a less good model for the Virgo IC LOSVD if the Virgo cluster is not fully virialized.

We then run a GMM analysis for each of these samples to isolate individual kinematic structures and trace the coexistence of the M87 halo/Virgo IC components. The GMM is a probabilistic model which assumes that a distribution of points, in our case the objects’ velocities, can be described as a linear combination of K independent Gaussian probability density functions (PDFs). The best set of Gaussians is selected by implementing the expectation-maximization (EM) algorithm, an iterative process that continuously updates the PDF parameters until convergence in the likelihood is reached. The process ends with the assignment of posterior probabilities, Γ, that each data point belongs to each of the K Gaussians. This analysis is based on the public scientific Python package scikit-learn (Pedregosa et al. 2011). However, our estimated Gaussian mixture distribution starts from a weighted sample of data, the weights, p<sub>gc</sub>, being the normalized probabilities to be a GC given its photometric properties (see Sec. 2.3.1). We then modified the python GMM routine accordingly, and the details on this will be given in Longobardi et al. (2018, in preparation). Here we only point out that the EM for mixture models consists of two steps. The first step, or E step, consists of calculating the expectation of the component assignments, Γ<sub>i,k</sub> for each data point given the model parameters. Thus, at the m<sup>th</sup> iteration of the EM process the posterior probabilities of the i<sup>th</sup> particle to belong to the k<sup>th</sup> Gaussian component become:

$$ \Gamma_{i,k}^{m} = \frac{p_k(x_i | \mu_k, \sigma_k) P_k^{m}}{\sum_{k=1}^{K} p_k(x_i | \mu_k, \sigma_k) P_k^{m}} \times P_{gc,i} $$

where the Gaussian density function (PDF) is written as $p(x_i) = \sum_{k=1}^{K} p_k(x_i | \mu_k, \sigma_k)P_k$, with $p_k(x | \mu_k, \sigma_k)$ being the individual mixture component centered on $\mu_k$, with a dispersion $\sigma_k$, and mixture weight $P_k$. The second step, or M step, consists of maximizing the expectations parameters, i.e., in updating the values $\mu_k, \sigma_k$, and $P_k$, as weighted parameters, the weights being $\Gamma_{i,k}$. In this way the EM process always takes into account the starting weight of each particle $p_{gc,i}$. Tests with the AIC show that the mixture that best describes the total GC LOSVD is a two component model in the first three elliptical bins, while in the outermost bin the GMM iden-
<figure><img src="image_6.png" /><figcaption>**Figure 6.** GMM analysis of the GC LOSVD carried out separately for the GCs located to the SE (Left Panels) and NW (Right Panels) of M87. From top to bottom the sample is divided in four elliptical annuli as given in Fig. 5, with N representing the total number of GCs used in the GMM fitting. A AIC analysis (2nd and 4th panels) prefers a two components Gaussian mixture out to R = 216 kpc, implying the contribution from the M87 galaxy halo and Virgo ICL to the total LOSVD at all distances out to this radius. For R > 216 kpc the best fit mixture is a model with one component tracing the Virgo cluster potential.</figcaption></figure>

tifies only one component. A comparison of the perfor-
mances of AIC and BIC scores was tested against mock
sets of data that mimic different LOSVDs and sample
sizes. AIC scores were found to be more reliable, hence
we only report these scores. Moreover, the stability of
the fitted parameters were tested with 1000 GMM runs
for different mock data sets and initialization values.

The results of our GMM analysis, the histograms of
the data together with the GMM best-fit model, and
relative AIC scores are shown in Fig. 6, while we sum-
marize the fitting Gaussian parameters in Table 1.

The GMM analysis shows that out to 216 kpc there
is the coexistence of a hotter (broader) and colder (nar-
rower) component, while outside this radius the nar-
rower Gaussian is not detectable, and the GC kinemat-
<table><caption>**Table 1.** Fitting parameters for the GMM best fit models for the M87 halo and Virgo ICL component and separately for the SE and NW regions</caption><thead><tr><th></th><th>SE<sub>M87,halo</sub></th><th>NW<sub>M87,halo</sub></th><th>SE<sub>ICL</sub></th><th>NW<sub>ICL</sub></th></tr></thead><tbody><tr><td colspan="5" style="text-align:center;">|R| ≤ 32 kpc</td></tr><tr><td>μ</td><td>1385 ± 30 km s<sup>-1</sup></td><td>1317 ± 22 km s<sup>-1</sup></td><td>1074 ± 403 km s<sup>-1</sup></td><td>1098 ± 300 km s<sup>-1</sup></td></tr><tr><td>σ</td><td>350 ± 21 km s<sup>-1</sup></td><td>307 ± 16 km s<sup>-1</sup></td><td>806 ± 285 km s<sup>-1</sup></td><td>424 ± 212 km s<sup>-1</sup></td></tr><tr><td>W</td><td>0.90 ± 0.07</td><td>0.90 ± 0.07</td><td>0.10 ± 0.04</td><td>0.10 ± 0.07</td></tr><tr><td colspan="5" style="text-align:center;">32 &lt; |R| ≤ 80 kpc</td></tr><tr><td>μ</td><td>1355 ± 34 km s<sup>-1</sup></td><td>1293 ± 31 km s<sup>-1</sup></td><td>1153 ± 421 km s<sup>-1</sup></td><td>865 ± 159 km s<sup>-1</sup></td></tr><tr><td>σ</td><td>371 ± 24 km s<sup>-1</sup></td><td>304 ± 22 km s<sup>-1</sup></td><td>729 ± 287 km s<sup>-1</sup></td><td>551 ± 112 km s<sup>-1</sup></td></tr><tr><td>W</td><td>0.88 ± 0.08</td><td>0.77 ± 0.7</td><td>0.12 ± 0.07</td><td>0.22 ± 0.06</td></tr><tr><td colspan="5" style="text-align:center;">80 &lt; |R| ≤ 216 kpc</td></tr><tr><td>μ</td><td>1319 ± 22 km s<sup>-1</sup></td><td>1282 ± 24 km s<sup>-1</sup></td><td>977 ± 182 km s<sup>-1</sup></td><td>986 ± 237 km s<sup>-1</sup></td></tr><tr><td>σ</td><td>234 ± 16 km s<sup>-1</sup></td><td>241 ± 17 km s<sup>-1</sup></td><td>680 ± 129 km s<sup>-1</sup></td><td>474 ± 167 km s<sup>-1</sup></td></tr><tr><td>W</td><td>0.82 ± 0.08</td><td>0.87 ± 0.09</td><td>0.18 ± 0.05</td><td>0.12 ± 0.06</td></tr><tr><td colspan="5" style="text-align:center;">|R| &gt; 216 kpc</td></tr><tr><td>μ</td><td>—</td><td>—</td><td>999 ± 135 km s<sup>-1</sup></td><td>833 ± 109 km s<sup>-1</sup></td></tr><tr><td>σ</td><td>—</td><td>—</td><td>634 ± 95 km s<sup>-1</sup></td><td>559 ± 77 km s<sup>-1</sup></td></tr><tr><td>W</td><td>—</td><td>—</td><td>1.0 ± 0.2</td><td>1.0 ± 0.2</td></tr></tbody></table>

ics are well described by one single, broad Gaussian. Within the uncertainties the narrow Gaussians describe the galaxy’s kinematics with a mean, systemic velocity consistent with $V_{LOS} = 1308$ km s$^{-1}$, and with a velocity dispersion that decreases as function of the distance, reaching a value of $\sigma_{LOS} \sim 230$ km s$^{-1}$ in the outermost bin. The broad Gaussians, instead, are consistent with the kinematics traced by the early-type dwarf population in the core of Virgo, the latter measured to peak at $V_{dE+dS0} = 1139 \pm 67$ km s$^{-1}$, with a velocity dispersion of $\sigma_{dE+dS0} = 649$ km s$^{-1}$ (Binggeli et al. 1993). The mean velocities and velocity dispersions of the ICGC component are consistent not only with those of the dwarf galaxies, but also with those of the more massive systems in the Virgo core (Binggeli et al. 1993). However, i) very broad, asymmetric wings, as the one we show in Fig. 5 and Fig. 6, are only observed for the LOSVD of the dwarf spheroidal galaxies (see Binggeli et al. 1993; Conselice et al. 2001), and ii) in the next Sections we will show that the photometric properties of the detected ICGCs are similar to what is measured for the dE system in Virgo, hence our comparison with the lower mass Virgo galaxies. The broad GC LOSVD component points toward a different kinematic behavior when measured separately in the two halves: SE of M87 the velocity dispersion is larger than in the NW region. However, more kinematic tracers will be needed to reduce the uncertainties associated with these results

(see Table 1). For clarity, the NW and SE GC velocity phase-spaces, together with the LOSVDs of the M87 halo and Virgo intra-cluster GCs are shown in Fig. 7.

Different LOSVDs for the galaxy halo and the ICL are predicted by cosmological analysis of structure formation. At the center of simulated clusters the total LOSVD splits in two Gaussians characterized by very different velocity dispersions (Dolag et al. 2010; Cui et al. 2014). One component is gravitationally bound to the galaxy and more spatially concentrated; the other is more diffuse and its high velocity dispersion reflects the satellites’ orbital distribution in the cluster gravitational potential. Using this framework, we can identify the narrow and the broad Gaussians traced by the GC kinematics with the M87 halo and Virgo ICL, the latter contributing ~13% of the GC sample within 216 kpc, beyond which the M87 halo component is no longer detected.

In Fig. 8 we plot again the comparison between the kinematics of the GC system and the Virgo dE population, however, this time, no separation between NW and SE fields is made. Interestingly, the kinematic behavior of the GC system exhibits a change in dynamics at a projected distance $R = 128$ kpc where the velocity dispersion reaches its minimum (black dots), and the mean velocity (green line) starts deviating from the systemic velocity of M87 (red line). Thus, interpreting such a distance as the truncation radius, $R_T$, of the
<figure><img src="image_7.png" /><figcaption>**Figure 7.** Top Panel: Projected velocity phase-space that compares the velocity distributions of the M87 halo GCs (green dots) and Virgo ICGCs (yellow dots) in the NW region. The histograms of the velocities, plotted on the side, compare the GC LOSVDs associated to the two dynamical components (green is M87 halo and yellow is ICL). Bottom Panel: Same as top panel, however the relation is shown for the SE region with respect to M87. The M87 systemic velocity ($V_{sys}$ = 1308 km s$^{-1}$) is plotted as red line. The vertical orange lines separate the phase-space in four elliptical bins, as given in the text.</figcaption></figure>

galaxy, would imply a total dark matter mass for M87, $M_{M87} = 4.6 \times 10^{13} M_{\odot}$. This is a rough estimate, that uses the relation:

$$R_T = R_{\text{Virgo}} \times \left( \frac{M_{\text{M87}}}{3 \times M_{\text{Virgo}}} \right)^{1/3}, \quad (1)$$

were $R_T = 127.5 \text{ kpc} \times ((1 - e)^{1/2})$ is the average distance within which systems are gravitationally bound to the galaxy, $R_{\text{Virgo}}$ is the distance between M87 and the center of the cluster (we have assumed $R_{\text{Virgo}} = 1^{\circ}$; Binggeli et al. 1993), and $M_{\text{Virgo}} = 4 \times 10^{14} M_{\odot}$ is the total Virgo dark matter mass (McLaughlin 1999). Moreover, we can also speculate that the different velocity dispersions measured for the IC component, higher in the south than in the north, are tracing the un-relaxed nature of the Virgo cluster. Young galaxy clusters still in the process of assembly are expected to show signs of dynamically non-relaxed substructures.

## 4. M87 HALO AND VIRGO INTRA-CLUSTER GC POPULATIONS

The statistical analysis we have carried out allows us to assign membership to either the M87 halo or the Virgo ICL. In this Section we further analyze the properties of the M87 galaxy and intra-cluster GCs that have spectroscopic information to investigate whether the different kinematic components also correspond to two different populations in terms of photometric properties and spatial distributions.

### 4.1. Photometric properties

Both the M87 halo and the Virgo intra-cluster GCs are characterized by a large spread in $(g' - i')_0$ colors, indicating that both dynamical components have two distinct, red (metal-rich) and blue (metal-poor), sub-populations of GCs. By fitting a double-Gaussian profile to the GC color distributions, we measure that the red and blue GC components cross at $(g' - i')_0 = 0.77 \text{ mag}$. The relative fraction of blue-to-red GCs differs, however:
<figure><img src="image_8.png" /><figcaption>**Figure 8.** Comparison between the LOS kinematics of the GCs in the NGVS pilot region and of the dwarf early-type galaxies within a 2° radius from the cluster's center. At a projected distance R = 128 kpc the GC velocity dispersion (black dots) reaches its minimum and the mean velocity (green line) starts deviating from the systemic velocity of M87 (red line). Beyond this radius, the GC kinematics traces the IC potential: the GC velocity dispersion rises till reaching the velocity dispersion values of the Virgo dE (purple crosses) at R ~ 320 kpc; at these distances the GC mean velocity merges into the mean velocity of the Virgo dwarfs (blue line). The mean velocity values are centered on V<sub>Sys</sub> = 1308 km s<sup>-1</sup>. The error bars show the uncertainties from Poissonian statistics.</figcaption></figure>

77% ± 10% of the Virgo ICGC population is blue, while
only 52% ± 3% of the GCs bound to the M87 halo are
metal-poor<fnref n="4" />. In Fig. 9, we show the $(g'-i')_0$ vs. $g'_0$ color-
magnitude diagram for both the M87 halo (green dots)
and the Virgo intra-cluster GCs (black crosses), sepa-
rately for the SE and NW fields, and as function of the
distance from the galaxy's center. The M87 halo sample
is compared in three different elliptical bins, while, given
the lower number of tracers, the ICGCs are analyzed
within two elliptical bins with R = 216 kpc representing
the boundary. Despite the incompleteness that charac-
terizes our sample of spectroscopically confirmed GCs
(higher for R > 216 kpc and with the NW side affected
the most), we see that the blue GC component of M87
is higher at larger distances, increasing from 40% close
to the center, to 70% in the outermost bin. This is not

true for the IC component that, independently from the
distance, is dominated by the blue GCs.

We also present the $(g' - r')_0$ vs. $(i' - z')_0$ color-color
diagrams for the M87 halo (Fig. 10) and the Virgo
intra-cluster GCs (Fig. 11). This analysis was moti-
vated by the recent finding of a correlation between
GC color-color relations and environment, and inter-
preted as driven by different stellar abundance ratios
(Powalka et al. 2016b, 2018, altough interpretations re-
lated to different underlying ages are also given, e.g.
Usher et al. (2015); Powalka et al. (2018) ). In this
study, the working sample was selected through a com-
bination of the $u^*i'K_s$ color properties and the compact-
ness of the objects, similar to what we have described in
Sec. 2.3.1. However, they only studied GCs whose mag-
nitudes were measured with an uncertainty $\sigma_{mag} < 0.06$
mag in all the bands, and also offset the $z'$ magnitudes
to account for known deviations of the SDSS-DR10 mag-
nitudes from the AB magnitude system (Powalka et al.
2016a,b). Hence, we consider the same constraints and
restrict our analysis to GCs with the same precision on
the magnitude values and apply the same shift to the $z'$
magnitudes. For each of the color-color diagrams we fit
and plot the relative Gaussian kernel density estimate
(kde) and locate its peak identifying the maximum val-
ues in $(g' - r')_0$ and $(i' - z')_0$ from the kde probability
density function.

The M87 halo GCs are shown in Fig. 10 (green dots), where we also flag the position of the blue and red class of GCs as identified based on their $(g' - i')_0$ values: green dots in red circles refer to the metal-rich population, while green dots in blue circles represent metal-poor GCs. Within $R \le 40$ kpc the SE and NW sides of M87 are characterized by GCs with similar colors, with the high-density peaks in the color-color space measured at $[(g' - r')_0, (i' - z')_0]_{SE} \sim [0.52, 0.15]$ $[0.51, 0.14]$, and $[(g' - r')_0, (i' - z')_0]_{NW} \sim [0.54, 0.14]$ $[0.52, 0.15]$. Similar are also the distributions of blue and red GCs in the color-color plane. At larger distances, i.e. $R > 40$ kpc the color-color diagram extends to bluer $(i' - z')_0$ colors, with a significant fraction of object with $(i' - z')_0 \le 0.1$. We note, however, that the average behavior differs between the SE and NW regions. NW of M87 the blue GCs define the high-density peak in the color-color diagram at $[(g' - r')_0, (i' - z')_0] \sim [0.50, 0.11]$. SE of M87 the blue GC population occupies a larger area and the high-density peak is measured at $[(g' - r')_0, (i' - z')_0] \sim [0.48, 0.14]$.

The same analysis is carried out for the ICGC sample,
(we use the same binning as in Fig.9). The results are
presented in Fig. 11, where this time black dots in red
and blue circles identify the metal-rich and metal-poor

<footnote n="4">⁴ Such a red-to-blue cluster ratios are consistent with those we would obtain integrating the extrapolated number density profiles for the red/blue populations related to the two dynamical components (see Sec.4.2)</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 9.** $(g' - i')_0$ vs $g'$ color-magnitude diagram for the M87 galaxy (green dots) and the Virgo intra-cluster (black crosses) GCs, at increasing distances from the galaxy's center (see legend) and for the SE and NW halves around M87. As a consequence of the smaller number of objects, the ICGC color-magnitude diagrams are compared in two different elliptical bins, with $R = 216$ kpc being the boarder. Both dynamical components have a blue and a red population of GCs. In the M87 halo the number of blue GCs increases as function of the distance: in the outermost bin it equals to ~70%.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 10.** $(g'-r')_0$ vs. $(i'-z')_0$ color-color diagrams and corresponding kde density contours for the GC populations associated to the M87 halo. **Left Panel:** The relation is shown for the GC population that lies in the SE side with respect to M87, and from top to bottom for GC sub-samples at increasing distances from the center of the galaxy. Independent of galactocentric distance, the high-density peaks are close in value as shown by the green lines. **Right Panel:** Same as in the left panel but for the sample of GCs in the NW half of M87. In this region the GCs show bluer $(i'-z')_0$ colors. In this diagram we flag the position of the blue and red class of GCs as identified based on their $(g' - i')_0$ values: green dots in red circles refer to the metal-rich population, while green dots in blue circles represent metal-poor GCs. The linear maximum-likelihood fits to the color-color relations, the $\Delta$ values, as well as the mean $(g' - z')_0$ colors for each GC sub-sample are given in the legend. The uncertainty on the magnitudes is less than 0.06 mag.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 11.** Same as in Fig. 10 but for the Virgo ICGCs. The ICGCs are characterized by bluer colors than the galaxy GCs (the green dashed line shows a representative peak at [(g' - r')_0, (i' - z')_0] ~ [0.51, 0.14]), with the bluest population found at R > 216 kpc. The magenta dashed lines show the color-color high-density peak for a sample of GCs associated to Virgo early-type dwarfs. The linear maximum-likelihood fits to the color-color relations, the Δ values, as well as the mean (g' - z')_0 colors for each GC sub-sample are given in the legend. The uncertainty on the magnitudes is less than 0.06 mag.</figcaption></figure>

populations of GCs, respectively. The ICGCs concentra-
tion in the color-color space differs from the one observed
for the M87 halo GCs: independently from the distance
a large fraction of ICGCs extend to (i' - z')_0 ≤ 0.1,
the distribution is dominated by the blue ICGCs, and
the few red GCs are bluer in the (g' - r')_0 vs (i' - z')_0
plane. This implies blue high-density peaks measured
at [(g' - r')_0, (i' - z')_0]SE ~ [0.44, 0.12] [0.48, 0.11] and
[(g' - r')_0, (i' - z')_0]NW ~ [0.47, 0.13] [0.49, 0.07]. These
are bluer than the one measured for the M87 halo GCs
of which we show a representative peak as green dashed
line. We also compared the ICGC color-color distribution
with the one measured for the 33 GCs that in Sec. 2.3.4
we identified as bound to Virgo early-type dwarfs. Their
peak (magenta dashed line in Fig. 11) is measured at

redder colors. Besides the offset in the color peaks, the
comparison between Fig. 10 and Fig.11 also shows dif-
ferent slopes for the M87 halo and Virgo intra-cluster
GCs color-color relation. A maximum-likelihood linear
fit to the (g' - r')_0 and (i' - z')_0 results into a galaxy
component that is usually steeper than the ICL. Also,
by computing the median deviation of the GC (i' - z')_0
colors from the fitted linear relation, Δ, the locus of the
M87 halo GCs is tighter (i.e., stronger color-color rela-
tion) than the one of the ICGCs (fitted parameters and
Δ values are given in the legend of Fig. 10 and Fig. 11).
Such properties can be interpreted in terms of different
metallicity distributions (Powalka et al. 2016b).

Finally, Peng et al. (2006) find that the mean (g' - z')_0
colors of the GCs in ACSVCS galaxies define a color-
magnitude sequence, such that more luminous/massive galaxies host, on average, redder GC systems. We than computed the mean $(g'-z')_0$ colors for the M87 halo and Virgo intra-cluster GCs and as function of the distance. The GC systems bound to the M87 halo show redder mean $(g'-z')_0$ values. Averaged across the distance, and with no distinction between NW and SE halves, we find $< (g'-z')_{0,M87} > = 1.07 \pm 0.01$, while for $R > 40$ kpc we find $< (g'-z')_0, M87 > = 1.00 \pm 0.01$. The Virgo ICGCs are constant, within the uncertainties, to a mean color $< (g'-z')_{0,IC} > 0.94 \pm 0.02$. Such a result is consistent with the galaxy accreted component dominating the outer regions and reinforce our conclusion that the ICL is built up by less massive systems than the one contributing to the size growth of the M87 halo.

## 4.2. Density profiles

We are now interested in studying the spatial distributions of the M87 halo and Virgo intra-cluster GCs to analyze their density profiles. This is because differences in their radial profiles may imply differences in the evolutionary paths of these two components as different progenitors are expected (Dolag et al. 2010; Cui et al. 2014) and observed (Longobardi et al. 2015a) to contribute differently to the resulting density distribution.

### 4.2.1. Spatial completeness

In order to analyze the GC density profiles we need to compute a completeness factor as function of the distance, $C_{spec}(R)$, that accounts for the spatial incompleteness of our spectroscopically selected sample of GCs. Such a completeness is estimated against the photometrically selected GCs, i.e. objects with high $p_{gc}$ probability (see Sec. 2.3.1), that we recall being 100% complete down to our limiting magnitude, $g' = 24.0$ mag. We then estimate the completeness function by computing the ratio between the number of photometric GCs with spectroscopic information, $N_{GC,spec}$, and the total number of GCs in the photometric sample, $N_{GC,phot}$, in each elliptical annulus<fnref n="5" />. Hence, the completeness factor as function of the distance from the M87’s center, $C_{spec}(R)$, can be written as:

$$C_{spec,j}(R) = \frac{N_{GC,spec,j}}{N_{GC,phot,j}},$$

where R is the average major axis distance of all GCs falling within each elliptical annulus. The subscript j in-

<figure><img src="image_12.png" /><figcaption>**Figure 12.** Spatial distribution of the velocity-selected sample of GCs (green dots are M87 halo GCs and yellow dots are ICGCs) superimposed on to the $g'$-band NGVS image of the central 2 by 2 square degree of Virgo. The NW field is more affected by spatial incompleteness. East is to the left, North is up.</figcaption></figure>

dicates the two different GC components, red and blue. Averaged over the elliptical bins the completeness factor amounts to $\langle C_{spec} \rangle = 0.1$, however the NW field is more affected by incompleteness than the SE region. This effect is measured at larger distances so that in the last bin we miss 40% more GCs NW to M87 as compared to the SE. This is due to observational strategies that observed the SE region more, but it is also related to the fact that the NW region has more photometric objects per unit area.

In Fig. 12 the sky positions of our sample of data (both M87 halo, green dots, and Virgo intra-cluster GCs, yellow dots) can be seen over-plotted on the NGVS $g'$-band image.

### 4.2.2. Density profile for the M87 halo and Virgo intra-cluster GCs

The density profiles for the M87 halo and intra-cluster GCs are constructed by binning our GC samples in elliptical bins, that have as major axis distances the same we used for our GMM analysis in Sec. 3 (see also orange lines in Fig. 5), and as position angle and ellipticity P.A. = -25.0°, and e = 0.4, respectively (Ferrarese et al. 2006; Janowiecki et al. 2010). Hence, for each of the two dynamical components of GCs and as function of

<footnote n="5">⁵ The completeness function $C_{spec}(R)$ is computed using the spectroscopic sample of objects retaining those GCs that we flagged as bound to other Virgo members. This is because the photometric sample of GCs we use to complete our observations does contain such a contribution, and in this way the ratio is still conserved.</footnote>
the distance, the density is computed as

$$ \rho_{\text{M87,j}}(R) = \frac{(N_{\text{obs,M87,j}}(R) - N_{\text{ov,j}}(R))}{A(R)} \times \frac{1}{C_{\text{spec,j}}(R)} \quad (2) $$

$$ \rho_{\text{ICL,j}}(R) = \frac{(N_{\text{obs,ICL,j}}(R) + N_{\text{ov,j}}(R))}{A(R)} \times \frac{1}{C_{\text{spec,j}}(R)} \quad (3) $$

where $N_{\text{obs,M87}}(R)/N_{\text{obs,ICL}}(R)$ is the observed number of M87 halo and Virgo intra-cluster GCs in the elliptical bin, respectively; $A(R)$ is the annulus area, estimated via Monte Carlo integration techniques if only a portion of the annulus intersects our FOV, and $C_{\text{spec}}(R)$ is the spectroscopic completeness factor. The observed numbers of clusters together with the completeness factor differ for the two, red and blue, populations of GCs, that in Eq. 2 and Eq. 3 are denoted with the subscript j. As can be seen in Fig. 6 the ICL and M87 halo velocity distributions overlap, and as result, the GMM algorithm assigns the ICGCs at low velocities relative to the galaxy systemic velocity to the galaxy halo component. We statistically quantified this effect comparing the M87 halo and IC velocity distributions in each elliptical bin. The LOSVDs of the two components were approximated by Gaussian functions with mean values and dispersions as given in Table 1. Hence, we calculated the fraction of ICGCs that lie inside the galaxy halo distribution as the area of overlap between the two curves. With this analysis we obtain a statistical estimate of the number of ICGCs, $N_{\text{ov}}(R)$, contained in the M87 halo velocity distribution. It is a function of R, however, averaged over the bins, 10% of the GC sample counted as part of the M87 halo is estimated to be associated with the ICL. For each of the radial bins, this contribution was subtracted from the M87 halo and added to the IC component. We emphasize that to build the density profiles associated with the red and blue GCs we consider the different contributions these populations have within the ICGC component and account for it when computing $N_{\text{ov}}$.

In Fig. 13 we show the comparison between the V-band surface brightness measured in the M87 region (black dots Janowiecki et al. 2010), the ICL-free M87 surface brightness (green dots), obtained subtracting from the measured $\mu_V$ the ICL surface brightness (see Sec. 4.2.3 for details about the adopted ICL surface brightness distribution), and the density profile of the M87 halo GCs (full dots) and Virgo ICGCs (stars). To this end we plot the logarithmic GC profile with its uncertainties defined as:

$$ \mu_{\text{GC}_1}(R) = -2.5 \log[\rho_1(R)] + k_{0t}, \quad (4) $$

$$ \sigma_{\mu_{\text{GC}_1}}(R) = \left| \frac{\partial \mu(R)}{\partial \rho_{\text{GC}_1}} \right| \times \sigma_{\rho_{\text{GC}_1}}(R), \quad (5) $$

where the subscript l denotes the two different GC components, and $k_{0t}$ is an arbitrary constant that we added for a better comparison with the light. The uncertainty on the density is given as $\sigma_{\rho_{\text{GC}_1}}(R) = 1/\sqrt{N_{\text{GC}_1}(R)} \times \rho_{\text{GC}_1}(R)$, with $N_{\text{GC}_1}(R)$ being the final number of GC counted in each radial bin. The density profiles are shown separately for the blue and red populations of GCs (blue and red filled dots for the M87 halo and blue and red stars for the Virgo intra-cluster GCs, respectively). Note that the normalization factors differ for the red and blue populations associated with both the M87 halo and ICL.

<figure><img src="image_13.png" /><figcaption>**Figure 13.** Top Panel: Azimuthally averaged density profiles as function of the major axis distance from M87 for the galaxy halo (circles) and IC (crosses) GCs, with the red and blue symbols identifying the red (metal-rich) and blue (metal-poor) populations of GCs. The M87 GCs show a distribution that is steeper than the IC counterpart, the latter being more extended, and shallower, however still centrally concentrated. The red GCs bound to the M87 halo (red circles) follow the ICL-free galaxy’s surface brightness profile (green dots). The blue GCs bound to the M87 halo (blue circles) are more extended deviating from the M87 light profile for R > 30 kpc. Bottom Panel: Difference in number of GCs between the NW and SE regions of M87 as function of the major axis distance from M87 (legend as in the top panel). A value of zero (dashed gray line) implies symmetry in number counts of GCs between the NW and SE fields. A significant deviation from such a symmetry is measured in the last bin for the ICGC population, implying a higher fraction of intra-cluster component in the NW region of M87.</figcaption></figure>

We find the red component of the M87 halo GCs agrees well with the galaxy’s light, following its surface brightness profile, while the blue component is more extended and flatter, in agreement with results from pre-
vious studies (e.g., Harris et al. 2009; Pota et al. 2013; Durrell et al. 2014). As discussed in the previous Section the total fraction of blue GCs in the M87 halo component amounts to 52% of the entire sample. However Fig. 13 shows that its relative contribution increases as function of the distance: close to the M87’s center, R ≤ 32 kpc, 40% of objects are blue, at intermediate distances, 32 < R ≤ 80 kpc, this fraction increases to 50%, and reaches 70% at larger distances. The ICGC density profile, which is flatter than the galaxy halo distribution, is still centrally concentrated. These findings are consistent with predictions from hydro-dynamical simulations, where the radial density profile of the component bound to the galaxy halo is much steeper than that of the diffuse IC component (Murante et al. 2004; Dolag et al. 2010; Cui et al. 2014).

Finally we compute the difference of GC counting for both the M87 halo and IC objects, and for the red and blue populations in the NW and SE fields with respect the M87’s center (Fig. 13), bottom panel). At all distances the galaxy halo component is consistent within the uncertainties with no asymmetry between the northern and southern fields (ΔN = 0, dashed line). The cluster component, instead, shows asymmetry in the number of ICGCs when the NW and SE fields are studied separately, implying a higher fraction of IC component in the NW field in the outermost bin. This observation is consistent with what is measured by Durrell et al. (2014) from the photometric sample and is now validated on the basis of the spectroscopically confirmed GCs.

### 4.2.3. The Virgo ICGCs $S_N$

The number of clusters per unit luminosity quantifies the GC specific frequency, $S_N$, that numerically is defined as

$$S_N = 2 \times N_{GC,TO} \times 10^{0.4(M_V+15)}, \quad (6)$$

for a system with absolute V-band magnitude $M_V$. $N_{GC,TO}$ is the number of GCs down to the GC luminosity function turnover magnitude, $\mu_{g,TO}$, and, multiplied by two, it is used to represent the total cluster population using a symmetric parametrization of the GC luminosity function (Harris & van den Bergh 1981). Assuming $\mu_{g,TO} = 24.0$ mag for the GC system at the center of Virgo (Jordán et al. 2007), $N_{ICGC,TO}$ in our surveyed area can be computed by fitting a power law to the measured ICGC density distribution and extrapolating inward such a profile. Doing so, we find that the expected number of ICGCs in the central $2^\circ \times 2^\circ$ of Virgo, and down to $g' \le 24.0$ mag, is $N_{ICGC,TO} = 2503 \pm 430$, where the error accounts for the uncertainty in the integrated density profile. The entire population of ICGCs can

then be retrieved multiplying by two $N_{ICGC,TO}$. However, such a choice is not free from caveats: we are assuming a constant GCLF turnover $g'$-band magnitude for the systems that contributed to the IC component that might result in an over-/underestimation of $N_{ICGC}$ for brighter and fainter $\mu_{g,TO}$, respectively. Thus, we computed the number of GCs we would add/lose if the 'true' GCLF would have magnitude cutoff in the range $\mu_{g,TO} \in [23.8, 24.2]$ with dispersion values $\sigma_{bright} = 1.3$ and $\sigma_{faint} = 0.8$. The latter represent the dispersion values found to be representative of the GCLF for bright and faint galaxies in Virgo, respectively. The range of $\mu_{g,TO}$ traces, instead, the uncertainty on the $g'$-band turnover magnitude for the GCLF of Virgo galaxies, estimated to be constant to a value $\mu_{g,TO} = 24.0 \pm 0.2$ (Jordán et al. 2007). Hence, we expect a total ICGC population of $N_{ICGC} = 5006 \pm 1021$, where the error now accounts for the 11% uncertainty due to the choice of the GCLF.

To estimate the total ICL luminosity we integrate the IC light profile traced by the ICPNs as given in Longobardi et al. (2018; in preparation) and in Longobardi et al. (2015a). Analyzing a sample of ICPNs in the inner 0.5 square degree of Virgo, these authors computed the V-band surface brightness profile for the IC component in units of $L_{\odot,V} \text{ pc}^{-2}$, as the scaled power-law fit to the IC PN surface density data (more details can be found in, Longobardi et al. 2015a). They found that $I(R)_{ICL} \propto R^{\gamma}$, with $\gamma = -0.79$ (Longobardi et al. 2018, in prep) or $\gamma \in [-0.34, -0.04]$ (Longobardi et al. 2015a). Integrating over the NGVS pilot region, this profile gives a total V-band luminosity $L_{ICL} = 4.2 \pm 0.1 \times 10^{10} L_{\odot}$, with the error on $L_{ICL}$ computed as the variation in luminosity values when the different light profiles are used. This leads to our estimate for the ICGC specific frequency that is $S_{N,ICL} = 10.2 \pm 4.8$. Peng et al. (2008) presented specific frequencies for the GC systems of 100 early-type galaxies in the ACS Virgo Cluster Survey and found that luminous early-type galaxies have higher $S_N$ than intermediate-luminosity early-type galaxies by a factor of ~2-3, with the latter class of objects ($-20.5 < M_V < -18$) generally having $S_N \sim 1.5$. Early-type dwarf galaxies have a large spread in $S_N$, characterized by an average value $S_N = 3.1$, however reaching estimates as high as $S_N = 14$. They also observed that nearly all dwarf galaxies with high GC fractions are within 1 Mpc of M87, hence the globular cluster fraction that we associate with the Virgo ICGC component is as high as the one measured for dwarf ellipticals that have yet to be subjected to tidal stripping or disruption (Peng et al. 2008).
5. DISCUSSION

5.1. *The Virgo intra-cluster population of GCs*

Durrell et al. (2014) presented a large-scale study of the distribution of GCs in Virgo, statistically accounting for background/foreground contamination to the GC sample. By analyzing the two-dimensional maps of the GC distribution they showed that the blue GCs represent an extended population, still significant beyond several effective radii from the hosting galaxy’s center. In the case of M87, they interpreted the GC population beyond 215 kpc as part of the intra-cluster component of Virgo, also suggested by a change in slope of the density profile of the blue GCs at major axis distances $R > 240$ kpc, and the evidence for a spatial asymmetry of GCs surrounding M87 for major axis distances larger than $20'$, with an excess of tracers in the NW region (mostly the blue population).

On the spectroscopic front the contamination by foreground MW stars has not allowed a clear identification of Virgo ICGCs. Strader et al. (2011) compiled a large catalog of GCs to study the halo of M87 and reported no evidence of an IC population of GCs neither in terms of a flattening of the GC density distribution nor in terms of a steep rise in the velocity dispersion. They find a low value of $\sigma \sim 300 \text{ km s}^{-1}$ in the galaxy’s outer regions with no difference in kinematics between the red and blue GCs. However, despite their efforts to reduce the MW star contamination the information they gathered was insufficient to disentangle contaminants and genuine GCs in the most ambiguous velocity range $V_{\text{LOS}} < 350 \text{ km s}^{-1}$, forcing them to classify all ambiguous objects as stars.

More recently, Ko et al. (2017) identified a sample of Virgo ICGCs that they associated with a main body of the cluster and with an infalling substructure. ICGCs with $V_{\text{LOS}} < \sim 400 \text{ km s}^{-1}$ were considered as representative of the infalling component, probably tracing systems that are dragged along by the M86 subgroup that is currently approaching the Virgo’s center with negative LOS velocities. We describe a similar phenomenon in Sec. 3.1, where the comparison between the NW and SE LOSVDs shows a higher fraction of dE galaxies with negative $V_{\text{LOS}}$ north of M87 (towards M86). However, our interpretation is that, in the process of falling, the M86 subgroup is contributing GCs to the IC component through the stripping of material coming from low mass galaxies that are being accreted into the Virgo core. These clusters hold the memory of the dynamics of their progenitors, however, they are not spatially associated with them, and we consider them as a genuine component of the ICL. We also note that the

ICGC LOSVD in Ko et al. (2017) does not extend below $V_{\text{LOS}} \sim -200 \text{ km s}^{-1}$, hence it is possible that they are missing some of the low velocity ICGCs due to their selection criteria, which were chosen to reduce foreground star contamination.

It is clear that the study of the Virgo ICGC component requires a proper accounting of contamination by foreground stars. In Sec. 2 we showed that the photometric information provided by the NGVS in the pilot region resulted in a clean photometric-based selection of GCs, thus, allowing us to avoid hard cuts in velocity and explore the entire Virgo cluster velocity range. We showed that the IC component, overlapping onto the one of the central galaxy, results in a total LOSVD with strong and quite asymmetric tails, also making the GC kinematics gradually deviate from the one of M87 at larger distances until it merges with the Virgo cluster values as traced by the early-type dwarf galaxies. The Virgo IC component, then, has a distinct dynamical populations of GCs representative of a dynamically young system. This is shown by the ICGC velocity distribution that does not show signs of point-symmetry—-independent of the distance, the ICGC velocity dispersion is higher in the SE region than in the NW region. Although more data is needed to make the latter statement statistically stronger, this result, together with the measured spatial asymmetry (see Sec.5.2) are consistent with the growing observational evidence of the ongoing assembly of Virgo. However, we note that Weil et al. (1997) argued that SE of M87 there is a spray of material as the result of a past interaction that could be responsible of the twisted and flaring outer isophotes to the SE of M87 as measured in the recent work by Mihos et al. (2017). Thus it is possible that there are multiple kinematic components in this region (see also Powalka et al. 2018), which manifest as a high velocity dispersion due to the limited number of tracers.

Finally, the Virgo ICGC component is not responsible for the high specific frequency measured for M87 ($\sim 3$ times higher than the average value observed for elliptical galaxies (Peng et al. 2008)). West et al. (1995), suggested that $\sim 70\%$ of the total M87 GC population was associated with cluster GCs. Although the results presented in this work do identify a Virgo ICGC population that overlaps with the M87 halo GCs, we find that the cluster only accounts for roughly 13% of the total number of GCs out to 216 kpc (see Table 1), a distance beyond which we do not predominantly trace the galaxy component anymore. Hence the presence of the Virgo ICGC population cannot explain the high M87 $S_N$.
## 5.2. Properties of the M87 halo and Virgo intra-cluster GCs and their connection to the systems formation

The kinematic separation between M87 halo and Virgo intra-cluster GCs allowed us in Sec. 4 to separately analyze their properties. We measured that both components have a red (metal-rich) and a blue (metal-poor) population of GCs, however while the majority of ICGCs are metal-poor, the red-to-blue GC ratio in M87 is close to 50%. On the other hand, the blue GC component of M87 increases as function of the radial distance till it dominates in the outermost bin. Different fractions of blue to red GCs also affect the $(g' - r')_0$ vs $(i' - z')_0$ color-color relation, that was recently shown to give insights in the evolution history of the host systems (Powalka et al. 2016b). We then find that the Virgo ICGCs peak on average at bluer colors than the M87 halo GCs. However the color-color regions occupied by the metal-poor GC component of the M87 halo population reach bluer colors moving towards large radii. In the NW region and for $R > 40$ kpc the peak for the M87 halo GCs reaches similar values to those representative of the IC population at those distances. A similar trend is observed for the ICGCs, that on the northern side and for $R > 216$ kpc are characterized by the bluest peak in the $(g' - r')_0$ vs $(i' - z')_0$ diagram.

We also have information on the density distribution of the M87 halo and intra-cluster GCs, with the latter measured to be shallower, although still centrally concentrated around M87. An excess of ICGCs in the NW region reveals the asymmetric nature of this component, while within the error the spatial distribution of the M87 halo does not deviate from point symmetry. Finally, we find a specific frequency for the Virgo ICGC population $S_{N,ICL} = 10.2 \pm 4.8$, assigning it to the category of high-$S_N$ systems.

The total number of GCs with respect to the parent stellar light, the relative fraction of blue-to-red objects, the GC color gradients as well as the different spatial distributions must relate to the dynamical evolution of the system that hosted these GCs. In what follows, we will review these properties and consider their implications for the formation and evolution of M87, the Virgo ICL and their GC populations.

**The Virgo cluster core**—Galaxy interactions, as well as tidal interactions between galaxies and the cluster potential are the main mechanisms expected to build the IC component over time. The evidence in Virgo points to an IC component made up of GCs that are stripped from their host systems orbiting in the cluster potential. 1) It is well established that blue GCs have a shallower

and more extended distribution than red GCs, and if so they are more likely to be stripped resulting in an IC component with a high fraction of blue to red GCs, and blue GC colors. Using N-body numerical simulations Ramos et al. (2015) analyzed the tidal removal of GCs in Virgo-like galaxies and found that on average halos lose ~16% and ~29% of their initial red and blue GC populations. On the other hand, observational works (e.g., Smith et al. 2013), also emphasize the impact of the orbital and spatial distributions on the stripping process. 2) The shallow and centrally concentrated density distribution that overlaps onto a steeper density is found to be representative of the simulated diffuse components at the center of clusters (Dolag et al. 2010; Cui et al. 2014) and representative of the low-mass progenitors that went to build up the ICL at low-redshift (Cooper et al. 2015). 3) Numerical simulations (Murante et al. 2004; Dolag et al. 2010) show that the stars that belong to the diffuse component did not lose memory of the kinematics of their progenitors, still reflecting the dynamics of the satellite. Thus, the kinematic result, revealing a GC system merging at large radii with the Virgo dwarf population, makes this class of galaxies the very likely progenitors of the IC component. This statement is also supported by the photometric properties of the Virgo ICGCs. The measured mean $(g' - z')_0$ colors are consistent with the one found to be representative of low mass systems (Peng et al. 2006). Furthermore, the $S_N$ value and the high blue-to-red ratio is consistent with the values measured for dwarf early-type galaxies found within 1 Mpc of the Virgo's center and mostly driven by the blue population of GCs (e.g., Peng et al. 2008; Mistani et al. 2016). The extreme environment these galaxies reside in favors a short star formation period (as indicated by their high-$[\alpha/Fe]$ abundances) making them very efficient at forming massive star clusters (Liu et al. 2016). Closer to the center (within 40 kpc from M87), however, dwarf galaxies have few or no GCs, suggesting that GCs have been stripped and are contributing to what today we see as the IC population. This picture is consistent with what is observed in deep optical surveys, where the color and morphology of the ICL stellar streams suggests that the objects getting shredded are low luminosity systems (e.g., Rudick et al. 2010; Mihos et al. 2017). However, it is important to note that we do observe a red GC fraction kinematically associated with the IC component, adding the result that among the satellites that in the past built up the diffuse light in Virgo there are also progenitors with $M_{prog} > 10^9 M_\odot$. This is because galaxies below such a mass have no measured red GC population (Peng et al. 2008).
In a recent work, Ferrarese et al. (2016) computed the expected total luminosity resulting from the disruption of galaxies in the Virgo core region, $L_{dis,V} = 7.7 \pm 2.0 \times 10^{10} L_{\odot}$. Hence, our estimate of the Virgo ICL luminosity implies that the IC component can contribute up to ~58% of the total ‘accreted luminosity’. Such an estimate assumes that the ICL only originates from galaxies that have been disrupted (and that had GCs), while it is likely that a fraction of the IC component has been stripped from still existing galaxies, or from galaxies that had no GCs.

Finally, as previously discussed, dE galaxies are the most likely progenitors of the Virgo ICL as traced by GCs. Thus, we can convert the total IC luminosity associated to the cluster component into a mass estimate, assuming as a mass-to-light-ratio for the ICL the one measured for these systems. Toloba et al. (2011) analyzed a sample of dE galaxies in Virgo and estimated an average dynamical mass-to-light-ratio in the I band, $\Upsilon_{I,dyn} = 2.0 \pm 0.04\Upsilon_{\odot}$. This implies a total mass for the IC component in the Virgo core $M_{ICC,tot} = 10.8 \pm 0.1 \times 10^{11} M_{\odot}$, that translates into a total stellar mass $M_{ICC,\star} = 6.6 \pm 5.0 \times 10^{11} M_{\odot}$ if $M_{tot}/M_{\star} = 1.6 \pm 1.2$ (Toloba et al. 2011)<fnref n="6" />. Needless to say, these results have important implications for the origin of the diffuse light in cluster cores, essential for a fair comparison with hydrodynamical cosmological simulations.

M87—Galaxies like M87 are thought to form via a two-phase formation scenario following which, to a first growth via rapid star formation fueled by the infall of cold gas or through major merger events, follows a series of minor merger episodes which enrich the galaxy outer regions of stars and make them grow efficiently in size. Numerical simulations predict that only very massive systems, like M87, are going to experience mergers up to ~5 times in their dynamical evolution history (De Lucia et al. 2006; De Lucia & Blaizot 2007), and their halos are built up through the accretion of both massive and small systems (Cooper et al. 2015), so that massive galaxies ($\log M_{\star} > 10.5 M_{\odot}$) are observed to assemble stellar mass predominantly in their halos (Eigenthaler et al. 2018).

The properties of the GCs associated with M87 likely reflect the evolutionary path taken by the galaxy itself. 1) The origin of the GC color bimodality was explained as the natural consequence of hierarchical merg-

ing. Following this interpretation, the metal-poor population of GCs represents the accreted component as a consequence of the accretion of smaller satellites, that adds on to the preexisting (in-situ) metal-rich GC population, the latter that formed along with the body of the galaxy itself (Côté et al. 1998; Kissler-Patig 2000; Mackey & Gilmore 2004). In this scenario, while the intrinsic (red) population of GCs dominates near the center of the galaxy, the ratio of metal-poor to metal-rich GCs increases with distance till the outer regions are characterized by an excess of blue clusters, as we measure for the M87 GCs. 2) On the other hand the higher fraction of red GCs measured moving closer to the center may also trace the accretion of more massive systems that due to dynamical friction would be dragged towards the central regions and there deposited their redder population of GCs.

Several questions remain unanswered. For example, can simulations of stripping of GCs from cluster dEs reproduce the different dynamical components of GCs that also differ in kinematics and density distributions? Would such a set of simulations be able to reproduce an IC component with a metal-rich population of clusters, and if so would they measure the same ratio of blue-to-red GCs we observe? What would be the dynamical/stellar mass associated with the diffuse component in simulated galaxy cluster cores? On the observational front, what causes the different peak in the $(g' - r')_0$ vs $(i' - z')_0$ color-color relation for the M87 halo GCs? Is it tracing the presence of the M87 crown NW of M87, that we know it is not a negligible perturbation of the galaxy halo properties, or is it tracing the accretion of different progenitors that differ in the SE and NW region? We could argue, indeed, that such a NW-SE division is the residual kinematical signature of progenitor galaxy groups having formed the center of Virgo, and we can speculate that if/when the M49 group will fall in, it will produce a similar structure for some Gyr in the N-S direction. Finally, it is possible that a fraction of the red GC component at the center of M87 can be traced back to the major merger event that caused the presence of the KDC (Emsellem et al. 2014). In this case, a different chemical composition could differentiate these objects from the red GC component that originated from other formation processes. Future projects focusing on GC stellar population properties derived through the analysis of the line-strengths of stellar absorption features in deep spectroscopy have the potential to untangle some of these issues.

## 6. SUMMARY

<footnote n="6">⁶ to transform V-band luminosities into I-band luminosities we have used $V - I = 1.03 \pm 0.04$ as calculated by van Zee et al. (2004) for dwarf galaxies.</footnote>
We have used the full photometric and spectroscopic information of the GCs in the NGVS Virgo core region (the 2° × 2° area around M87) to identify and separate the M87 halo GC population from the Virgo ICGCs, and study their properties separately. This study can be summarized as follows:

* We implemented the extreme deconvolution (XD) algorithm to model as a sum of Gaussians the multi-dimensional color-concentration-magnitude space occupied by the objects in the NGVS, and reduce the contamination by foreground stars. This allowed us to explore the entire Virgo velocity range, usually affected by hard cuts in velocity for $V_{LOS} < 500$ km s⁻¹ (Fig. 2 and Fig. 3).

* The GC kinematics reveal a system that, as function of the distance from the M87's center, deviates from the galaxy's kinematics. On both sides of M87, the GC system shows a decrease in the mean velocity and an increase in the velocity dispersion until it merges (for distances larger than 320 kpc) with the Virgo cluster potential, as traced by the dwarf early-type galaxies in Virgo (Fig. 5).

* Based on a Gaussian Mixture Model (GMM) analysis, the mixture that best describes the total LOSVD is a two component model, that we then interpret as arising from the Virgo's IC and M87's halo components. M87 halo and ICL overlap out to 216 kpc. Further out, the GMM identifies only one component with mean velocity and velocity dispersion consistent with tracing the ICL (Fig. 6). This analysis, run separately for the sample of data in the SE and NW regions of M87, shows an IC component that is not characterized by point symmetry: the SE region shows higher values of velocity dispersion with respect to the NW counterpart. This finding can be interpreted as the result of an IC component that is not relaxed. However, due to the limited number of GC velocity measurements at these extreme distances, the level of significance is within $1\sigma$.

* The analysis of the $(g'-i')_0$ vs. $g'$ color-magnitude diagram shows that both M87 galaxy and Virgo ICGC populations have a blue (metal-poor, $(g'-i')_0 < 0.77$) and a red (metal-rich, $(g'-i')_0 \ge 0.77$) component of GCs. However the IC population is richer in blue GCs, the latter representing ~70% of the IC sample. The fraction of blue GCs in the M87 component is ~50% of the sample (Fig. 9). On the other hand, the blue GC component of

the M87 halo increases as function of the distances reaching 70% in the outermost bin.

* The ICGCs are characterized by bluer colors than the galaxy halo GC population as shown by the study of the $(g'-r')_0$ vs. $(i'-z')_0$ color-color diagrams. Within the M87 halo population the GCs in the NW region and for major axis distances $40\,\text{kpc} < R \le 216\,\text{kpc}$ are bluer, and they may be driven by the presence of a recent accretion event that in that region we know is contributing 60% of the M87 halo light (Longobardi et al. 2015b). Within the IC population, the GCs in the NW region and for major axis distances $R > 216\,\text{kpc}$ are the bluest (Fig. 10 and Fig. 11).

* The intra-cluster GC population has a shallower density profile than the M87 galaxy GC population, consistent with results from theoretical simulations of diffuse stellar component (ICL) around central cluster galaxies (Dolag et al. 2010). The different slopes can also be interpreted as due to the different origins of the two components: steeper profiles come from the accretion of more massive systems at higher redshifts, while the shallower and more extended ICL component would be the result of less massive systems whose accretion happened at lower redshifts (e.g., Cooper et al. 2015). This is consistent with the measured mean $(g'-z')_0$ colors that are consistent with the one found to be representative of low luminosity/mass systems (Peng et al. 2006). We find evidence for an asymmetry in the number of ICGCs when the NW and SE fields are studied separately, implying a higher fraction of IC component in the NW field for distances larger than $R > 216\,\text{kpc}$ (Fig. 13).

* The population of ICGCs falls under the category of high-$S_N$ systems with a measured specific frequency $S_{N,ICL} = 10.2 \pm 4.8$. This value is consistent with the one observed for the dwarf early-type galaxies that reside within 1 Mpc from the cluster center, supporting a stripping / disruption scenario for the ICGCs.

* The total V-band luminosity relative to the ICL at the center of Virgo is $L_{ICL} = 4.2 \pm 0.1 \times 10^{10} L_\odot$. This value translates into a total (dark plus luminous) and stellar masses $M_{ICC,tot} = 10.8 \pm 0.1 \times 10^{11} M_\odot$, and $M_{ICC,\star} = 6.6 \pm 5.0 \times 10^{11} M_\odot$, respectively.

These results show that the GCs successfully trace the galaxy halo-ICL interface, and their properties allow
us to understand the formation and evolution of these
two components. This work depicts a picture in which
the M87 halo and Virgo intra-cluster GCs have different
progenitors, with the latter most likely representing the
Virgo accreted component from low and intermediate
mass dwarf-ellipticals.

ACKNOWLEDGMENTS

AL is thankful to E. Emsellem and L. Sales for con-
structive and useful discussion on the analysis and in-
terpretation of the data. EWP and SL acknowledge
support from the National Natural Science Foundation
of China through Grant No. 11573002. EWP also ac-
knowledges travel support from the Chinese Academy
of Sciences South America Center for Astronomy (CAS-
SACA) that was essential for developing the photomet-
ric classification techniques used in this paper. SE and
AJ acknowledge support from project IC120009 “Mil-
lennium Institute of Astrophysics (MAS)” of the Mil-
lennium Science Initiative, Chilean Ministry of Econ-
omy. AJ acknowledges additional support from project
BASAL CATA PFB-06. THP. acknowledges support
through the FONDECYT Regular Project No. 1161817
and the BASAL Center for Astrophysics and Associ-
ated Technologies (PFB-06). THP and AL acknowl-
edge ECOS-Sud/CONICYT project C15U02. H.-X.Z.
acknowledges a support from the CAS Pioneer Hundred
Talents Program. C.L. acknowledges the NSFC grants
11673017, 11203017, 11433002. C.L. is supported by
Key Laboratory for Particle Physics, Astrophysics and

Cosmology, Ministry of Education. This paper is based
on observations obtained with MegaPrime/MegaCam, a
joint project of CFHT and CEA/IRFU, at the Canada-
France-Hawaii Telescope (CFHT) which is operated by
the National Research Council (NRC) of Canada, the
Institut National des Sciences de l'Univers of the Cen-
tre National de la Recherche Scientifique (CNRS) of
France, and the University of Hawaii. This research used
the facilities of the Canadian Astronomy Data Centre
operated by the National Research Council of Canada
with the support of the Canadian Space Agency. Ob-
servations reported here were obtained at the MMT
Observatory, a joint facility of the University of Ari-
zona and the Smithsonian Institution. MMT telescope
time was granted, in part, by NOAO, through the Tele-
scope System Instrumentation Program (TSIP). TSIP
is funded by NSF. Data presented in this paper were
obtained at the Anglo-Australian Telescope, which is
operated by the Australian Astronomical Observatory.
Data presented herein were obtained at the W.M. Keck
Observatory, which is operated as a scientific partner-
ship among the California Institute of Technology, the
University of California and the National Aeronautics
and Space Administration. The Observatory was made
possible by the generous financial support of the W.M.
Keck Foundation. This research made use of Astropy
(http://www.astropy.org), a community-developed core
Python package for Astronomy (Astropy Collaboration,
2018).

REFERENCES

Agnello, A., Evans, N. W., Romanowsky, A. J., & Brodie, J. P. 2014, MNRAS, 442, 3299

Aguerri, J. A. L., Gerhard, O. E., Arnaboldi, M., et al. 2005, AJ, 129, 2585

Alam, S., Albareti, F. D., Allende Prieto, C., et al. 2015, ApJS, 219, 12

Alamo-Martínez, K. A., & Blakeslee, J. P. 2017, ApJ, 849, 6

Arnaboldi, M., Freeman, K. C., Mendez, R. H., et al. 1996, ApJ, 472, 145

Arnaboldi, M., Gerhard, O., Aguerri, J. A. L., et al. 2004, ApJL, 614, L33

Arp, H. C. 1967, Astrophys. Lett., 1, 1

Ashman, K. M., & Zepf, S. E. 1992, ApJ, 384, 50

Ashman, K. M., & Bird, C. M. 1993, AJ, 106, 2281

Barbary, K., Aldering, G., Amanullah, R., et al. 2012, ApJ, 745, 31

Barbosa, C. E., Arnaboldi, M., Coccato, L., et al. 2018, A&A, 609, A78

Battaglia, G., Helmi, A., Morrison, H., et al. 2005, MNRAS, 364, 433

Baum, W. A. 1955, PASP, 67, 328

Beasley, M. A., Baugh, C. M., Forbes, D. A., Sharples, R. M., & Frenk, C. S. 2002, MNRAS, 333, 383

Bender, R., Saglia, R. P., & Gerhard, O. E. 1994, MNRAS, 269, 785

Bender, R., Kormendy, J., Cornell, M. E., & Fisher, D. B. 2015, ApJ, 807, 56

Bekki, K., & Yahagi, H. 2006, MNRAS, 372, 1019

Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393

Bertin, E., Mellier, Y., Radovich, M., et al. 2002, Astronomical Data Analysis Software and Systems XI, 281, 228

Binggeli, B., Sandage, A., & Tammann, G. A. 1985, AJ, 90, 1681

Binggeli, B., Tammann, G. A., & Sandage, A. 1987, AJ, 94, 251
Binggeli, B., Popescu, C. C., & Tammann, G. A. 1993, A&AS, 98, 275

Blakeslee, J. P., Jordán, A., Mei, S., et al. 2009, ApJ, 694, 556

Blakeslee, J. P., Cho, H., Peng, E. W., et al. 2012, ApJ, 746, 88

Boselli, A., Boissier, S., Cortese, L., & Gavazzi, G. 2008, ApJ, 674, 742-767

Boselli, A., Voyer, E., Boissier, S., et al. 2014, A&A, 570, A69

Bovy, J., Hennawi, J. F., Hogg, D. W., et al. 2011, ApJ, 729, 141

Brodie, J. P., & Huchra, J. P. 1991, ApJ, 379, 157

Brodie, J. P., & Strader, J. 2006, ARA&A, 44, 193

Brodie, J. P., Romanowsky, A. J., Strader, J., et al. 2014, ApJ, 796, 52

Cantiello, M., & Blakeslee, J. P. 2007, ApJ, 669, 982

Carter, D., & Dixon, K. L. 1978, AJ, 83, 574

Coccato, L., Arnaboldi, M., & Gerhard, O. 2013, MNRAS, 436, 1322

Conselice, C. J., Gallagher, III, J. S., & Wyse, R. F. G. 2001, ApJ, 559, 791

Contini, E., De Lucia, G., Villalobos, Á., & Borgani, S. 2014, MNRAS, 437, 3787

Cooper, A. P., Gao, L., Guo, Q., et al. 2015, MNRAS, 451, 2703

Côté, P., Marzke, R. O., & West, M. J. 1998, ApJ, 501, 554

Côté, P., Marzke, R. O., West, M. J., & Minniti, D. 2000, ApJ, 533, 869

Côté, P., McLaughlin, D. E., Hanes, D. A., et al. 2001, ApJ, 559, 828

Côté, P., West, M. J., & Marzke, R. O. 2002, ApJ, 567, 853

Côté, P., McLaughlin, D. E., Cohen, J. G., & Blakeslee, J. P. 2003, ApJ, 591, 850

Côté, P., Blakeslee, J. P., Ferrarese, L., et al. 2004, ApJS, 153, 223

Cui, W., Murante, G., Monaco, P., et al. 2014, MNRAS, 437, 816

De Lucia, G., Springel, V., White, S. D. M., Croton, D., & Kauffmann, G. 2006, MNRAS, 366, 499

De Lucia, G. & Blaizot, J. 2007, MNRAS, 375, 2

Dilday, B., Bassett, B., Becker, A., et al. 2010, ApJ, 715, 1021

Dolag, K., Murante, G., & Borgani, S. 2010, MNRAS, 405, 1544

Drinkwater, M. J., Jones, J. B., Gregg, M. D., & Phillipps, S. 2000, PASA, 17, 227

Durrell, P. R., Ciardullo, R., Feldmeier, J. J., Jacoby, G. H., & Sigurdsson, S. 2002, ApJ, 570, 119

Durrell, P. R., Côté, P., Peng, E. W., et al. 2014, ApJ, 794, 103

Emsellem, E., Krajnović, D., & Sarzi, M. 2014, MNRAS, 445, L79

Eigenthaler, P., Puzia, T. H., Taylor, M. A., et al. 2018, ApJ, 855, 142

Feldmeier, J. J., Ciardullo, R., Jacoby, G. H., & Durrell, P. R. 2004a, ApJ, 615, 196

Ferguson, A. M. N., Gallagher, J. S., & Wyse, R. F. G. 1998, AJ, 116, 673

Ferrarese, L., Côté, P., Jordán, A., et al. 2006, ApJS, 164, 334

Ferrarese, L., Côté, P., Cuillandre, J.-C., et al. 2012, ApJS, 200, 4

Ferrarese, L., Côté, P., Sánchez-Janssen, R., et al. 2016, ApJ, 824, 10

Firth, P., Drinkwater, M. J., & Karick, A. M. 2008, MNRAS, 389, 1539

Forbes, D. A., Brodie, J. P., & Grillmair, C. J. 1997, AJ, 113, 1652

Forbes, D. A., & Forte, J. C. 2001, MNRAS, 322, 257

Forbes, D. A. 2017, MNRAS, 472, L104

Ford, H. C., & Butcher, H. 1979, ApJS, 41, 147

Forte, J. C., Martinez, R. E., & Muzzio, J. C. 1982, AJ, 87, 1465

Gal-Yam, A., Maoz, D., Guhathakurta, P., & Filippenko, A. V. 2003, AJ, 125, 1087

Gavazzi, G., Boselli, A., Vílchez, J. M., Iglesias-Paramo, J., & Bonfanti, C. 2000, A&A, 361, 1

Gavazzi, G., Boselli, A., Donati, A., Franzetti, P., & Scodeggio, M. 2003, A&A, 400, 451

Gebhardt, K., & Kissler-Patig, M. 1999, AJ, 118, 1526

Georgiev, I. Y., Puzia, T. H., Goudfrooij, P., & Hilker, M. 2010, MNRAS, 406, 1967

Gerhard, O. E. 1993, MNRAS, 265, 213

Gerhard, O., Arnaboldi, M., Freeman, K. C., et al. 2005, ApJL, 621, L93

Gnedin, O. Y. 2003, ApJ, 589, 752

Gonzalez, A. H., Zabludoff, A. I., & Zaritsky, D. 2005, ApJ, 618, 195

Gonzalez, A. H., Zaritsky, D., & Zabludoff, A. I. 2007, ApJ, 666, 147

Grossauer, J., Taylor, J. E., Ferrarese, L., et al. 2015, ApJ, 807, 88

Hanes, D. A., Côté, P., Bridges, T. J., et al. 2001, ApJ, 559, 812

Harris, W. E., & van den Bergh, S. 1981, AJ, 86, 1627

Harris, W. E. 1991, ARA&A, 29, 543

Harris, W. E., Kavelaars, J. J., Hanes, D. A., Pritchet, C. J., & Baum, W. A. 2009, AJ, 137, 3314
Harris, W. E., Harris, G. L. H., & Alessi, M. 2013, ApJ, 772, 82

Harris, W. E., Ciccone, S. M., Eadie, G. M., et al. 2017, ApJ, 835, 101

Haşegan, M., Jordán, A., Côté, P., et al. 2005, ApJ, 627, 203

Hartke, J., Arnaboldi, M., Longobardi, A., et al. 2017, A&A, 603, A104

Hartke, J., Arnaboldi, M., Gerhard, O., et al. 2018, arXiv:1805.03092

Hilker, M., Infante, L., Vieira, G., Kissler-Patig, M., & Richtler, T. 1999, A&AS, 134, 75

Hudson, M. J., & Robison, B. 2018, MNRAS

Janowiecki, S., Mihos, J. C., Harding, P., et al. 2010, ApJ, 715, 972

Jordán, A., Côté, P., West, M. J., & Marzke, R. O. 2002, ApJL, 576, L113

Jordán, A., West, M. J., Côté, P., & Marzke, R. O. 2003, AJ, 125, 1642

Jordán, A., Côté, P., West, M. J., et al. 2004, AJ, 127, 24

Jordán, A., McLaughlin, D. E., Côté, P., et al. 2007, ApJS, 171, 101

Kartha, S. S., Forbes, D. A., Spitler, L. R., et al. 2014, MNRAS, 437, 273

Kelson, D. D., Zabludoff, A. I., Williams, K. A., et al. 2002, ApJ, 576, 720

Kim, S., Rey, S.-C., Jerjen, H., et al. 2014, ApJS, 215, 22

King, I. R. 1966, AJ, 71, 64

Kissler-Patig, M. 2000, Reviews in Modern Astronomy, 13, 13

Ko, Y., Hwang, H. S., Lee, M. G., et al. 2017, ApJ, 835, 212

Kormendy, J., Fisher, D. B., Cornell, M. E., & Bender, R. 2009, ApJS, 182, 216

Kruijssen, J. M. D. 2014, Classical and Quantum Gravity, 31, 244006

Kundu, A., & Whitmore, B. C. 2001, AJ, 121, 2950

Laporte, C. F. P., White, S. D. M., Naab, T., & Gao, L. 2013, MNRAS, 435, 901

Larsen, S. S., Brodie, J. P., Huchra, J. P., Forbes, D. A., & Grillmair, C. J. 2001, AJ, 121, 2974

Leaman, R., VandenBerg, D. A., & Mendel, J. T. 2013, MNRAS, 436, 122

Lee, M. G., Park, H. S., & Hwang, H. S. 2010, Science, 328, 334

Lim, S., Peng, E. W., Duc, P.-A., et al. 2017, ApJ, 835, 123

Liu, Y., Zhou, X., Ma, J., et al. 2005, AJ, 129, 2628

Liu, C., Peng, E. W., Côté, P., et al. 2015, ApJ, 812, 34

Liu, C., Peng, E. W., Toloba, E., et al. 2015, ApJL, 812, L2

Liu, Y., Peng, E. W., Blakeslee, J., et al. 2016, ApJ, 818, 179

Longobardi, A., Arnaboldi, M., Gerhard, O., et al. 2013, A&A, 558, A42

Longobardi, A., Arnaboldi, M., Gerhard, O., & Hanuschik, R. 2015, A&A, 579, A135

Longobardi, A., Arnaboldi, M., Gerhard, O., & Mihos, J. C. 2015, A&A, 579, L3

Mackey, A. D., & Gilmore, G. F. 2004, MNRAS, 355, 504

McDonald, M., Courteau, S., Tully, R. B., & Roediger, J. 2011, MNRAS, 414, 2055

McLaughlin, D. E. 1999, ApJL, 512, L9

Mei, S., Blakeslee, J. P., Côté, P., et al. 2007, ApJ, 655, 144

McLaughlin, D. E., Harris, W. E., & Hanes, D. A. 1994, ApJ, 422, 486

Merritt, D. 1984, ApJ, 276, 26

Mihos, J. C., Harding, P., Feldmeier, J., & Morrison, H. 2005, ApJL, 631, L41

Mihos, J. C., Harding, P., Feldmeier, J. J., et al. 2017, ApJ, 834, 16

Misgeld, I., & Hilker, M. 2011, MNRAS, 414, 3699

Mistani, P. A., Sales, L. V., Pillepich, A., et al. 2016, MNRAS, 455, 2323

Montes, M., Trujillo, I., Prieto, M. A., & Acosta-Pulido, J. A. 2014, MNRAS, 439, 990

Moore, B., Lake, G., & Katz, N. 1998, ApJ, 495, 139

Mould, J. R., Oke, J. B., & Nemec, J. M. 1987, AJ, 93, 53

Muñoz, R. P., Puzia, T. H., Lançon, A., et al. 2014, ApJS, 210, 4

Murante, G., Arnaboldi, M., Gerhard, O., et al. 2004, ApJL, 607, L83

Naab, T., Johansson, P. H., & Ostriker, J. P. 2009, ApJL, 699, L178

Nulsen, P. E. J. & Bohringer, H. 1995, MNRAS, 274, 1093

Oldham, L. J., & Evans, N. W. 2016, MNRAS, 462, 298

Oser, L., Ostriker, J. P., Naab, T., Johansson, P. H., & Burkert, A. 2010, ApJ, 725, 2312

Ohyama, Y., & Hota, A. 2013, ApJL, 767, L29

Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, Journal of Machine Learning Research, 12, 2825

Peng, E. W., Ford, H. C., & Freeman, K. C. 2004, ApJ, 602, 705

Peng, E. W., Jordán, A., Côté, P., et al. 2006, ApJ, 639, 95

Peng, E. W., Jordán, A., Côté, P., et al. 2008, ApJ, 681, 197-224

Peng, E. W., Jordán, A., Blakeslee, J. P., et al. 2009, ApJ, 703, 42

Peng, E. W., Ferguson, H. C., Goudfrooij, P., et al. 2011, ApJ, 730, 23

Peng, E. W., & Lim, S. 2016, ApJL, 822, L31

Perrett, K. M., Stiff, D. A., Hanes, D. A., & Bridges, T. J. 2003, ApJ, 589, 790
Powalka, M., Lançon, A., Puzia, T. H., et al. 2016, ApJS, 227, 12

Powalka, M., Puzia, T. H., Lançon, A., et al. 2016, ApJL, 829, L5

Powalka, M., Lançon, A., Puzia, T. H., et al. 2017, ApJ, 844, 104

Powalka, M., Puzia, T. H., Lançon, A., et al. 2018, ApJ, 856, 84

Pota, V., Graham, A. W., Forbes, D. A., et al. 2013, MNRAS, 433, 235

Puchwein, E., Springel, V., Sijacki, D., & Dolag, K. 2010, MNRAS, 406, 936

Puzia, T. H., Kissler-Patig, M., Brodie, J. P., & Huchra, J. P. 1999, AJ, 118, 2734

Puzia, T. H., Kissler-Patig, M., Thomas, D., et al. 2004, A&A, 415, 123

Puzia, T. H., Perrett, K. M., & Bridges, T. J. 2005, A&A, 434, 909

Puzia, T. H., Kissler-Patig, M., Thomas, D., et al. 2005, A&A, 439, 997

Roediger, J. C., Ferrarese, L., Côté, P., et al. 2017, ApJ, 836, 120

Sand, D. J., Graham, M. L., Bildfell, C., et al. 2011, ApJ, 729, 142

Schuberth, Y., Richtler, T., Hilker, M., et al. 2010, A&A, 513, A52

Smith, R., Sánchez-Janssen, R., Fellhauer, M., et al. 2013, MNRAS, 429, 1066

Sparks, W. B., Ford, H. C., & Kinney, A. L. 1993, ApJ, 413, 531

Stetson, P. B., Vandenberg, D. A., Bolte, M., Hesser, J. E., & Smith, G. H. 1989, AJ, 97, 1360

Strader, J., Brodie, J. P., Cenarro, A. J., Beasley, M. A., & Forbes, D. A. 2005, AJ, 130, 1315

Ramos, F., Coenda, V., Muriel, H., & Abadi, M. 2015, ApJ, 806, 242

Richtler, T. 2013, 370 Years of Astronomy in Utrecht, 470, 327

Romanowsky, A. J., Strader, J., Brodie, J. P., et al. 2012, ApJ, 748, 29

Rudick, C. S., Mihos, J. C., & McBride, C. 2006, ApJ, 648, 936

Rudick, C. S., Mihos, J. C., Frey, L. H., & McBride, C. K. 2009, ApJ, 699, 1518

Rudick, C. S., Mihos, J. C., Harding, P., et al. 2010, ApJ, 720, 569

Sánchez-Janssen, R., Ferrarese, L., MacArthur, L. A., et al. 2016, ApJ, 820, 69

Santos, M. R. 2003, Extragalactic Globular Cluster Systems, 348

Schroder, L. L., Brodie, J. P., Kissler-Patig, M., Huchra, J. P., & Phillips, A. C. 2002, AJ, 123, 2473

Schuberth, Y., Richtler, T., Bassino, L., & Hilker, M. 2008, A&A, 477, L9

Seigar, M. S., Graham, A. W., & Jerjen, H. 2007, MNRAS, 378, 1575

Spengler, C., Côté, P., Roediger, J., et al. 2017, ApJ, 849, 55

Strader, J., Romanowsky, A. J., Brodie, J. P., et al. 2011, ApJS, 197, 33

Taylor, M. A., Puzia, T. H., Harris, G. L., et al. 2010, ApJ, 712, 1191

Toloba, E., Guhathakurta, P., Peletier, R. F., et al. 2014, ApJS, 215, 17

Toloba, E., Boselli, A., Cenarro, A. J., et al. 2011, A&A, 526, A114

Toloba, E., Li, B., Guhathakurta, P., et al. 2016, ApJ, 822, 51

Thomas, D., Maraston, C., Bender, R., & Mendes de Oliveira, C. 2005, ApJ, 621, 673

Toledo, I., Melnick, J., Selman, F., et al. 2011, MNRAS, 414, 602

Tonini, C. 2013, ApJ, 762, 39

Tully, R. B., & Shaya, E. J. 1984, ApJ, 281, 31

Usher, C., Forbes, D. A., Brodie, J. P., et al. 2015, MNRAS, 446, 369

van Zee, L., Barton, E. J., & Skillman, E. D. 2004, AJ, 128, 2797

Veale, M., Ma, C.-P., Greene, J. E., et al. 2017, MNRAS, 471, 1428

Ventimiglia, G., Arnaboldi, M., & Gerhard, O. 2011, A&A, 528, A24+

Vilchez-Gómez, R. 1999, in Astronomical Society of the Pacific Conference Series, Vol. 170, The Low Surface Brightness Universe, ed. J. I. Davies, C. Impey, & S. Phillips, 349

Watson, D. F. & Conroy, C. 2013, ApJ, 772, 139

Weil, M. L., Bland-Hawthorn, J., & Malin, D. F. 1997, ApJ, 490, 664

West, M. J., Cote, P., Jones, C., Forman, W., & Marzke, R. O. 1995, ApJL, 453, L77

West, M. J., Jordán, A., Blakeslee, J. P., et al. 2011, A&A, 528, A115

Williams, B. F., Ciardullo, R., Durrell, P. R., et al. 2007, ApJ, 654, 835

Yoshida, M., Yagi, M., Okamura, S., et al. 2002, ApJ, 567, 118

Zepf, S. E., Beasley, M. A., Bridges, T. J., et al. 2000, AJ, 120, 2928
Zhang, H.-X., Peng, E. W., Côté, P., et al. 2015, ApJ, 802, 30

Zhu, L., Long, R. J., Mao, S., et al. 2014, ApJ, 792, 59

Zhang, H.-X., Puzia, T. H., Peng, E. W., et al. 2018, ApJ, 858, 37

Zibetti, S., White, S. D. M., Schneider, D. P., & Brinkmann, J. 2005, MNRAS, 358, 949

Zwicky, F. 1951, PASP, 63, 61

Zwicky, F. 1952, PASP, 64, 242