# The Next Generation Virgo Cluster Survey. XXXIII. Stellar Population Gradients in the Virgo Cluster Core Globular Cluster System

Youkyung Ko, Eric W. Peng, Patrick Côté, Laura Ferrarese, Chengze Liu, Alessia Longobardi, Ariane Lançon, Roberto P. Muñoz, Thomas H. Puzia, Karla A. Alamo-Martínez, et al.

► To cite this version:

Youkyung Ko, Eric W. Peng, Patrick Côté, Laura Ferrarese, Chengze Liu, et al. The Next Generation Virgo Cluster Survey. XXXIII. Stellar Population Gradients in the Virgo Cluster Core Globular Cluster System. Astrophys.J., 2022, 931 (2), pp.120. 10.3847/1538-4357/ac63cf. hal-03659540

HAL Id: hal-03659540

https://hal.science/hal-03659540

Submitted on 17 Mar 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# The Next Generation Virgo Cluster Survey. XXXIII. Stellar Population Gradients in the Virgo Cluster Core Globular Cluster System

Youkyung Ko<fnref n="1" />,<fnref n="2" />,<fnref n="3" /> ID, Eric W. Peng<fnref n="2" />,<fnref n="3" /> ID, Patrick Côté<fnref n="4" /> ID, Laura Ferrarese<fnref n="4" /> ID, Chengze Liu<fnref n="5" /> ID, Alessia Longobardi<fnref n="6" />, Ariane Lançon<fnref n="7" /> ID, Roberto P. Muñoz<fnref n="8" /> ID, Thomas H. Puzia<fnref n="8" /> ID, Karla A. Alamo-Martínez<fnref n="9" /> ID, Laura V. Sales<fnref n="10" /> ID, Felipe Ramos-Almendares<fnref n="11" />,<fnref n="12" />,<fnref n="13" /> ID, Mario G. Abadi<fnref n="11" />,<fnref n="14" /> ID, Myung Gyoon Lee<fnref n="15" /> ID, Ho Seong Hwang<fnref n="1" />,<fnref n="15" />,<fnref n="16" /> ID, Nelson Caldwell<fnref n="17" /> ID, John P. Blakeslee<fnref n="18" /> ID, Alessandro Boselli<fnref n="19" /> ID, Jean-Charles Cuillandre<fnref n="20" /> ID, Pierre-Alain Duc<fnref n="7" /> ID, Susana Eyheramendy<fnref n="21" />,<fnref n="22" /> ID, Puragra Guhathakurta<fnref n="23" /> ID, Stephen Gwyn<fnref n="4" /> ID, Andrés Jordán<fnref n="21" />,<fnref n="22" /> ID, Sungsoon Lim<fnref n="24" /> ID, Rubén Sánchez-Janssen<fnref n="25" /> ID, and Elisa Toloba<fnref n="26" /> ID

## Abstract

We present a study of the stellar populations of globular clusters (GCs) in the Virgo Cluster core with a homogeneous spectroscopic catalog of 692 GCs within a major-axis distance $R_{maj} = 840$ kpc from M87. We investigate radial and azimuthal variations in the mean age, total metallicity, [Fe/H], and $\alpha$-element abundance of blue (metal-poor) and red (metal-rich) GCs using their co-added spectra. We find that the blue GCs have a steep radial gradient in [Z/H] within $R_{maj} = 165$ kpc, with roughly equal contributions from [Fe/H] and [ $\alpha$/Fe], and flat gradients beyond. By contrast, the red GCs show a much shallower gradient in [Z/H], which is entirely driven by [Fe/H]. We use GC-tagged Illustris simulations to demonstrate an accretion scenario where more massive satellites (with more metal- and $\alpha$-rich GCs) sink further into the central galaxy than less massive ones, and where the gradient flattening occurs because of the low GC occupation fraction of low-mass dwarfs disrupted at larger distances. The dense environment around M87 may also cause the steep [ $\alpha$/Fe] gradient of the blue GCs, mirroring what is seen in the dwarf galaxy population. The progenitors of red GCs have a narrower mass range than those of blue GCs, which makes their gradients shallower. We also explore spatial inhomogeneity in GC abundances, finding that the red GCs to the northwest of M87 are slightly more metal-rich. Future observations of GC stellar population gradients will be useful diagnostics of halo merger histories.

*Unified Astronomy Thesaurus concepts: Globular star clusters (656); Chemical abundances (224); Virgo Cluster (1772); Stellar populations (1622); Giant elliptical galaxies (651); Galaxy clusters (584)*

## 1. Introduction

The assembly of the most massive galaxies in the universe is a story of continual merging and accretion of gas, stars, and dark matter over a Hubble time. These behemoths have long ago stopped forming significant numbers of their own stars

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ Korea Astronomy and Space Science Institute, 776 Daedeok-daero, Yuseong-Gu, Daejeon 34055, Republic of Korea; ykko@kasi.re.kr, peng@pku.edu.cn</footnote>

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

<footnote n="3">³ Kavli Institute for Astronomy and Astrophysics, Peking University, Beijing 100871, People's Republic of China</footnote>

<footnote n="4">⁴ Herzberg Astronomy and Astrophysics Research Centre, National Research Council of Canada, 5071 W. Saanich Road, Victoria, BC V9E 2E7, Canada</footnote>

<footnote n="5">⁵ Department of Astronomy, School of Physics and Astronomy, and Shanghai Key Laboratory for Particle Physics and Cosmology, Shanghai Jiao Tong University, Shanghai 200240, People's Republic of China</footnote>

<footnote n="6">⁶ Dipartimento di Fisica G. Occhialini, Università degli Studi di Milano Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy</footnote>

<footnote n="7">⁷ Université de Strasbourg, CNRS, Observatoire astronomique de Strasbourg, UMR7550, F-67000 Strasbourg, France</footnote>

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

<footnote n="9">⁹ Department of Astronomy, DCNE, Universidad de Guanajuato, Callejon de Jalisco s/n, 36023, Guanajuato, Mexico</footnote>

<footnote n="10">¹⁰ University of California Riverside, 900 University Ave., Riverside, CA 92521, USA</footnote>

<footnote n="11">¹¹ CONICET—Universidad Nacional de Córdoba, Instituto de Astronomía Teórica y Experimental (IATE), Laprida 854, X5000BGR, Córdoba, Argentina</footnote>

<footnote n="12">¹² Instituto de Investigación Multidisciplinar en Ciencia y Tecnología, Universidad de La Serena, Raúl Bitrán 1305, La Serena, Chile</footnote>

<footnote n="13">¹³ Departamento de Física y Astronomía, Universidad de La Serena, Av. Cisternas 1200 N, La Serena, Chile</footnote>

<footnote n="14">¹⁴ Observatorio Astronómico, Universidad Nacional de Córdoba, Laprida 854, X5000BGR, Córdoba, Argentina</footnote>

<footnote n="15">¹⁵ Astronomy Program, Department of Physics and Astronomy, Seoul National University, Gwanak-gu, Seoul 151-742, Republic of Korea</footnote>

<footnote n="16">¹⁶ SNU Astronomy Research Center, Seoul National University, 1 Gwanak-ro, Gwanak-gu, Seoul 08826, Republic of Korea</footnote>

<footnote n="17">¹⁷ Center for Astrophysics, Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA</footnote>

<footnote n="18">¹⁸ NSF's NOIRLab, 950 N. Cherry Avenue, Tucson, AZ 85719, USA</footnote>

<footnote n="19">¹⁹ Aix Marseille Univ, CNRS, CNES, LAM, Marseille, France</footnote>

<footnote n="20">²⁰ AIM, CEA, CNRS, Université Paris-Saclay, Université Paris Diderot, Sorbonne Paris Cité, Observatoire de Paris, PSL University, F-91191 Gif-sur-Yvette Cedex, France</footnote>

<footnote n="21">²¹ Facultad de Ingeniería y Ciencias, Universidad Adolfo Ibáñez, Av. Diagonal las Torres 2640, Peñalolén, Santiago, Chile</footnote>

<footnote n="22">²² Millennium Institute of Astrophysics, Santiago, Chile</footnote>

<footnote n="23">²³ Department of Astronomy and Astrophysics, University of California, Santa Cruz, UCO/Lick Observatory, 1156 High Street, Santa Cruz, CA 95064, USA</footnote>

<footnote n="24">²⁴ Department of Astronomy, Yonsei University, 50 Yonsei-ro, Seodaemun-gu, Seoul 03722, Republic of Korea</footnote>

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

<footnote n="26">²⁶ Department of Physics, University of the Pacific, 3601 Pacific Avenue, Stockton, CA 95211, USA</footnote>

<footnote>Received 2021 November 2; revised 2022 March 31; accepted 2022 April 1; published 2022 June 1</footnote>
in situ, but they are still adding stars to their halos through the accretion of lower-mass galaxies (Oser et al. 2010; Johansson et al. 2012). When massive galaxies are situated at the centers of galaxy clusters, their stellar halos ("cD envelopes") join seamlessly into the larger accreted population sometimes referred to as "intracluster light" (Murante et al. 2004, 2007; Purcell et al. 2007; Contini et al. 2014). While the relatively low surface brightnesses of these halo stars make them difficult to study, they contain a wealth of information on the progenitor components that assembled to form the parent galaxy. In the most massive galaxies, the accreted halo stars may actually be the dominant stellar population (Purcell et al. 2007; Oser et al. 2010; Cooper et al. 2013).

The characterization of extragalactic halo stellar populations is made challenging by their low surface brightnesses, usually well below that of the sky. Deep imaging surveys of massive galaxies have revealed extended stellar halos in cluster galaxies (Bernstein et al. 1995; Gonzalez et al. 2000; Feldmeier et al. 2004; Mihos et al. 2005, 2017; Krick & Bernstein 2007; Rudick et al. 2010), some showing significant fine structure expected from merging (Duc et al. 2015; Bílek et al. 2020), and some have been resolved into individual stars (Mouhcine et al. 2005a, 2005b, 2007; Williams et al. 2012). Obtaining further information on the nature of these stellar halos is more difficult, as even colors are subject to uncertainties with sky subtraction (e.g., Mihos et al. 2013). Recent spectroscopic surveys have been able to derive the integrated ages, metallicities, and α-element abundances of the brighter regions ($R < 3$ effective radii) of the stellar halos of some massive early-type galaxies (e.g., Spolaor et al. 2010; Greene et al. 2015; Gu et al. 2018, 2018; Iodice et al. 2019; Fensch et al. 2020).

Numerous studies have used discrete tracers like planetary nebulae (PNe; e.g., Hui et al. 1993; Arnaboldi et al. 1996, 1998; Ciardullo et al. 2004; Buzzoni et al. 2006; Ciardullo 2010; Cortesi et al. 2013; Longobardi et al. 2013, 2015a) and globular clusters (GCs; e.g., Beasley et al. 2000; Peng et al. 2006; Hwang et al. 2008; Liu et al. 2011; Pastorello et al. 2015; Usher et al. 2015; Lim et al. 2017; Longobardi et al. 2018b; Ko et al. 2018, 2019, 2020; Powalka et al. 2018) with success to elucidate the stellar content of halo populations. In particular, colors of old GCs are indicative of metallicities, which have been broadly used to trace the mass assembly history of host galaxies. The GCs in massive galaxies commonly show bimodal color distributions, with the blue and red GCs corresponding to metal-poor and metal-rich subpopulations, respectively. The properties of red GCs show tighter correlations with the underlying stars in their host galaxies in terms of their spatial distributions (e.g., Park & Lee 2013; Wang et al. 2013) and kinematics (e.g., Peng et al. 2004; Pota et al. 2013). Several studies have suggested a scenario where the red (metal-rich) GCs were formed with the bulk of the stars during dissipative mergers, while the blue (metal-poor) GCs were formed in small galaxies and have been continuously accreted onto massive galaxies (e.g., Côté et al. 1998; Lee et al. 2010b).

We can better understand these possibilities by combining chemical and spatial information. Metallicity gradients, a common feature of galaxies, inform us about the role of the gravitational potential in the regulation of star formation and chemical enrichment, but they also highlight the role of major mergers in erasing these gradients (e.g., Spolaor et al. 2010). Gas-poor minor merging can also create metallicity gradients,

as higher-mass galaxies with more metal-rich stars have a tendency to sink further into the host galaxy owing to dynamical friction (Amorisco 2019). For example, in the Milky Way's stellar halo, the existence or not of abundance gradients has been important in shaping our view of how the halo may have formed. The lack of a metallicity gradient for Milky Way GCs beyond 8 kpc led Searle & Zinn (1978) to propose that the outer halo was formed from the accretion of low-mass fragments.

The color gradients of GC systems have been explored in a relatively large sample of early-type galaxies in the Virgo and Fornax Clusters (Liu et al. 2011), finding that both the blue (metal-poor) and red (metal-rich) GC systems exhibit radial color gradients, showing equal strengths within the uncertainties. This result holds even when gradient measurements are extended to $R \sim 5R_e-8R_e$, as in the literature compilation of Forbes & Remus (2018), but flattening of the gradients of metal-poor GCs was found at $R \gtrsim 8R_e$.

Another chemical parameter for understanding the formation history of galaxies is the ratio between α-element abundances and iron, [α/Fe]. This provides information about the star formation timescale, i.e., a high [α/Fe] value indicates rapid star formation and early quenching. As an example, the correlation of [α/Fe] in Milky Way stars with their locations in the Galaxy has been used to understand the formation of thin and thick disks (e.g., Recio-Blanco et al. 2014). For GCs in massive early-type galaxies, it is known that they have a wide range of [α/Fe] values and often have supersolar [α/Fe] ratios (e.g., Puzia et al. 2005; Park et al. 2012), which indicates that they were formed over a short timescale. Studies of the spatial variation of [α/Fe] in extragalactic GCs, however, are rare, but they have the potential to reveal in detail how galaxy halos are assembled.

Recently, there have been attempts to understand the characteristics of the GC populations in the context of cosmological simulations. The E-MOSAICS project (Pfeffer et al. 2018; Kruijssen et al. 2019; Reina-Campos et al. 2019) combined the MOSAICS semianalytic models for star cluster formation and evolution with the EAGLE hydrodynamical simulation for galaxy formation, including 25 Milky Way-like disk galaxies. They studied the formation and coevolution of galaxies and their star clusters in terms of chemical abundances, progenitors, and accretion epochs of the GCs. In another approach, Ramos-Almendares et al. (2020) investigated the GC populations in nine simulated Virgo-like galaxy clusters. They tagged the GCs to the simulated galaxies at their infall time and traced the kinematics and spatial distributions of the GCs. Other groups have also started to tag GCs in cosmological simulations to study the spatially resolved properties of GC systems (Trujillo-Gomez et al. 2021; Chen & Gnedin 2022). These developments have made it possible to use the spatial gradients and distributions of GCs in meaningful comparisons to simulations.

The Virgo Cluster is an excellent target to investigate the mass assembly history of the central galaxy because of its proximity. It is a dynamically young system that shows many substructures of member galaxies both spatially and kinematically (Binggeli et al. 1985, 1987). In particular, the core region of the Virgo Cluster shows diffuse stellar features (Mihos et al. 2005, 2017) and extended structures of ionized gas filaments (Kenney et al. 2008; Boselli et al. 2018) and H I gas (Yoshida et al. 2002; Oosterloo & van Gorkom 2005)
around galaxies, indicative of strong galaxy interactions. The most massive galaxy at the center of the Virgo Cluster, M87, has a weak cD envelope, and kinematic studies of GCs and PNe have been done to understand the connection between its faint stellar halo and the intracluster light (Longobardi et al. 2015a; Ko et al. 2017; Longobardi et al. 2018a, 2018b). In addition, the GCs and PNe in M87 also constitute kinematic substructures related to recent mergers (Strader et al. 2011; Romanowsky et al. 2012; Longobardi et al. 2015b; Oldham & Evans 2016; Lambert et al. 2020). Therefore, M87 and the Virgo Cluster are a good laboratory to study the ongoing assembly of the most massive galaxy in a cluster core.

In addition to the many aforementioned GC kinematic studies, the chemistry and ages of GCs have also been studied from photometry and spectroscopy in Virgo ellipticals (e.g., Cohen et al. 1998, 2003). Powalka et al. (2018) reported a substructure of relatively young GCs located ~11.5 south of M87, possibly related to past mergers or accretions. More recently, Villaume et al. (2020) extended the study of chemical abundances in GCs out to ~140 kpc from M87, showing the promise of stellar population studies in extragalactic GCs across a range of radius. They do not find any significant iron abundance gradients of both metal-poor and metal-rich populations within R < 40 kpc and conclude that the flat metallicity gradient of metal-poor GCs could be a sign that some metal-poor GCs were formed in situ at high redshift. Moreover, they found that the metal-poor GC population has a higher α-element abundance than the metal-rich GC population in both the inner (R < 40 kpc) and outer (R > 40 kpc) halos, which is not expected by current cosmological simulations.

In this study, we investigate mean trends of the stellar populations of the GCs in the core region of the Virgo Cluster out to ~840 kpc (2°9), as a function of major-axis distance from M87, using a new homogeneous spectroscopic catalog. This survey mainly covers a 3° × 2° field including the massive elliptical galaxies M87, M86, and M84. This paper is organized as follows. The GC catalog used in this study and the sample selection are described in Section 2. In Section 3, we explain how we measure ages and metallicities of the GC subpopulations. We present the main results about the age and metallicity distributions of the GCs in Section 4. We discuss the assembly history of M87 and the Virgo Cluster core in Section 5. In Section 6, we summarize the results. We adopt a distance to the Virgo Cluster of 16.5 Mpc (Mei et al. 2007; Blakeslee et al. 2009), so 1' corresponds to 4.91 kpc.

## 2. Data and GC Sample Selection

### 2.1. Sample Definition

We produced a new homogeneous catalog of the GCs in the Virgo core region by combining two spectroscopic data sets. The first data set was part of a spectroscopic follow-up program to the Next Generation Virgo Cluster Survey (NGVS; Ferrarese et al. 2012). The NGVS is an extensive deep imaging survey using the MegaCam imager on the Canada-France-Hawaii Telescope (CFHT). It images a 104 deg² field of the Virgo Cluster with five optical filters $u^*g'r'i'z'$. The images were taken with excellent seeing conditions (<0.6 for i-band imaging), and many GCs in the Virgo Cluster are marginally resolved in the NGVS images. The i'-band “fuzziness” index ($\Delta i_{4-8}$), defined as the difference between 4-pixel and 8-pixel i'-band aperture-corrected magnitude, has larger values for the

GCs than for foreground stars on average (Durrell et al. 2014). In addition, a $K_s$-band imaging survey was carried out with WIRCam on the CFHT (NGVS-IR; Muñoz et al. 2014) for a 2° × 2° region around M87. Muñoz et al. (2014) showed that the $u^*i'K_s$ color-color diagram allows us to clearly distinguish Virgo GCs from contaminants such as foreground stars and background galaxies. We identified GCs in the Virgo core based on a “fuzziness” index and various $u^*g'r'i'z'K_s$ color combinations of point sources, adopting the “extreme deconvolution” (xD) classification algorithm (Bovy et al. 2011). GCs were selected using xD, but also using a broader set of color and morphology cuts, and targeted for spectroscopy. This program obtained spectra of 776 GCs in the Virgo Cluster core, some of which have been used in previous studies (Zhu et al. 2014; Li et al. 2020).

Second, we used spectra taken from Ko et al. (2017), who presented a spectroscopic catalog of the GCs in a central 2° × 2° field of the Virgo Cluster, covering the vicinity of M87, M84, and M86. They obtained 910 optical spectra of GC candidates. We reclassified their sample using the $u^*i'K_s$ color-color diagram for the objects in the NGVS-IR field and confirmed 212 GCs in their sample. There are 88 sources in common between the NGVS follow-up program and Ko et al. (2017). Finally, we combined these two samples to construct a combined spectroscopic catalog of 903 GCs in total in the Virgo core.

Both programs were taken with the MMT/Hectospec (Fabricant et al. 2005) with a 270 mm⁻¹ grating having a dispersion of 1.2 Å pixel⁻¹. Data were obtained in the wavelength range 3650–9200 Å. Data reduction was done with the SAO pipeline, following the description in Caldwell et al. (2011), and flux calibration was done following the steps in Fabricant et al. (2008). We used the xcsao task in the IRAF RVSAO package (Kurtz & Mink 1998) to measure heliocentric radial velocities of targets, taking 10 different templates including stars, GCs, and galaxies. Combining radial velocity, color, and size information, we classify the objects into four groups: GCs, UCDs, foreground stars, and background galaxies. The details of the data reduction, flux calibration, heliocentric radial velocity measurements, and target classification will be described in a forthcoming paper about the spectroscopic catalog of the GCs (Y. Ko et al., 2022, in preparation). The median signal-to-noise ratio (S/N) of the spectra of GC candidates at 5000 Å is S/N ~9, where the Poisson noise, sky noise, and readout noise are considered to calculate spectral noise as a function of wavelength.

We select 726 out of the 903 total GCs, excluding the GCs that belong to other Virgo galaxies in the core. The remaining GCs are either GCs that belong to M87 or intracluster GCs (IGCs), a population governed by the galaxy cluster gravitational potential (Williams et al. 2007; Firth et al. 2008; Lee et al. 2010a; Durrell et al. 2014). We did this because our interest for this paper is mainly to investigate the assembly history of M87 in the cluster environment. There are 154 Virgo galaxies that have radial velocity information in the survey region. We used both spatial and kinematic criteria to separate the GCs bound to Virgo galaxies, following Longobardi et al. (2018b): (1) $R < 10R_e$, where $R_e$ is the effective radius of a given galaxy; and (2) $|v_r - v_{gal}| < 3\sigma_{gal}$ where $v_r$, $v_{gal}$, and $\sigma_{gal}$ are the radial velocities of GCs, galaxy systemic velocity, and galaxy central velocity dispersion, respectively. The radial velocities of galaxies are compiled from various literature
sources (NASA/IPAC Extragalactic Database (NED) 2019; Ahn et al. 2014; Kim et al. 2014; Ferrarese et al. 2020). The velocity dispersions are from the GOLDMine database (Gavazzi et al. 2003) compiling data from McElroy (1995) and Scodeggio et al. (1998). For the galaxies where velocity dispersions are not known, mostly faint dwarf galaxies, we adopted $\sigma_{gal} = 50$ km s⁻¹.

We expect that these cuts are very conservative and remove from our sample nearly all the GCs that belong to other galaxies. As an example, for larger galaxies, the effective radius of the GC system ($R_{e,GC}$) may be as big as ~5 times that of the stars (e.g., Kartha et al. 2016; Caso et al. 2019), so our cut to eliminate GCs within $10R_{e,stars}$ corresponds to a cut of $2R_{e,GC}$, which results in ~85% of GCs being removed, leaving ~15% to be potential contaminants. The velocity cut, however, removes all GCs within $\pm 3 \times$ the velocity dispersion from the mean galaxy velocity. If the GC line-of-sight velocity distribution (LOSVD) is roughly Gaussian, then this selection removes all but 0.27% of the GCs. Combined with the spatial cut, only 0.04% of the GC system could potentially contaminate our sample. Our spectroscopic limit ($g' < 22.5$ mag) and a Gaussian GC luminosity function shape ($\langle g \rangle \approx 24$ mag, $\sigma \sim 1.3$ mag; Jordán et al. 2007) mean that we only observe the brightest ~12% of the GCs in any massive galaxy. Therefore, the fraction of GCs that could enter our sample is roughly $5 \times 10^{-5}$. The massive nearby galaxies M86 and M84 have $\approx 3000$ GCs each (Peng et al. 2008). The analysis above therefore leads us to expect that only 0.3 GCs from these galaxies could contaminate our sample. Of course, this analysis assumes that the GC spatial distribution follows a Sérsic profile at large radii and that the LOSVD is Gaussian in the wings. These simple assumptions may not hold in practice, but any GC that is significantly detached (spatially and kinematically) from its host is probably on its way to being accreted by the M87 system and arguably should be included in our analysis.

Figure 1(a) shows the radial velocities of 726 GCs as a function of major-axis distance from M87. The major-axis distance $R_{maj}$ is calculated as follows:

$$R_{maj} = [(x/(1 - e))^2 + y^2]^{1/2}, \quad (1)$$

where $x = \Delta R.A. \cos(PA) - \Delta decl. \sin(PA)$

and $y = \Delta R.A. \sin(PA) + \Delta decl. \cos(PA)$,

with a position angle (PA) of 155° and an ellipticity $e$ of 0.4 for M87 (Ferrarese et al. 2006; Janowiecki et al. 2010). Most GCs have radial velocities consistent with a mean value of 1255 km s⁻¹, similar to the systemic velocity of M87 (1284 km s⁻¹; Cappellari et al. 2011). In addition, we identify another group of GCs at $30' < R_{maj} < 150'$, that appears to be a broad feature in the velocity-position phase-space diagram connecting M87 with M86. The mean radial velocity of these GCs is $v_r = -66$ km s⁻¹, which is about 1370 km s⁻¹ lower than the systemic velocity of M87.

Figure 1(b) shows the spatial distribution of the GCs in our study. The GCs are strongly concentrated around M87, although some of them are uniformly distributed in the survey region beyond $R_{maj} = 65'$. The GCs with low radial velocities are mostly located in the region west of M87 (between it and M86). Note that the low-velocity GCs that are strongly concentrated on M86 (systemic velocity of -224 km s⁻¹) are considered to belong to it and are not plotted in this figure. For this analysis, we only used the 692 GCs that have higher radial

<figure><img src="image_4.png" /><figcaption>Figure 1. (a) Radial velocity vs. major-axis distance from M87 for the M87 GCs and IGCs. The black circles and blue crosses represent the GCs in the Virgo main body and the low-velocity substructure of the GCs, respectively. The red plus signs indicate the locations of M86 and M84. The systemic velocities of M87 (1284 km s⁻¹), M86 (-224 km s⁻¹), and M84 (1017 km s⁻¹) are marked with horizontal dashed lines (Cappellari et al. 2011). The black solid line indicates the criterion for separating out the low-velocity objects. (b) Spatial distribution of the M87 GCs and IGCs. The center of M87 is at (0, 0), and the symbols are the same as in panel (a). The dotted outlines show the survey region of this study.</figcaption></figure>

velocities in order to focus on the main structure of the Virgo Cluster, assuming that the low-velocity GCs constitute an infalling group associated with M86.

## 2.2. GC Subsamples

Many studies have confirmed the existence of color bimodality for GCs in M87 (Whitmore et al. 1995; Elson & Santiago 1996; Kundu et al. 1999; Larsen et al. 2001; Peng et al. 2006, 2009; Tamura et al. 2006; Harris 2009), a feature typical of GCs in massive early-type galaxies. Since we expect that the blue and red GC subpopulations have different formation origins, we divided our GC sample into blue and red, using their $(g' - i')_0$ colors, to search for any differences in their stellar populations.

We performed Gaussian mixture modeling (Muratov & Gnedin 2010) on the $(g' - i')_0$ color distribution of the GCs in a $2^\circ \times 2^\circ$ region around M87, identified in the NGVS photometric survey (see Section 2.1) within radial bins $5'$ in width, assuming heteroscedastic bimodal distributions. Figure 2(a) shows the $(g' - i')_0$ color of the photometric sample as a function of
<figure><img src="image_5.png" /><figcaption>**Figure 2.** (a) $(g' - i')_0$ color vs. major-axis distance from M87 for the GC candidates within a 2° × 2° region around M87, identified in the NGVS photometric survey. The blue squares and red circles represent the mean $(g' - i')_0$ colors of the blue and red GC candidates within each radial bin, respectively, estimated from the GMM. The error bars indicate the Gaussian widths of the two components. The green solid lines indicate the color criterion for dividing the GCs into blue and red GCs. (b) Same as panel (a), but for the GCs without the low-velocity GCs, spectroscopically confirmed in this study.</figcaption></figure>

major-axis distance from M87. We adopted the color at which GCs have a 50% probability of belonging to either component as the $(g' - i')_0$ color that divides the two components. This color criterion changes as a function of major-axis distance. We adopted a broken linear fit for the dividing color for GCs with $R_{maj} < 20'$ as follows:

$$ (g' - i')_0 = 0.8490 - 0.0760 \times R_{maj}, \quad (2a) $$

$$ (g' - i')_0 = 0.8048 - 0.0013 \times R_{maj}, \quad (2b) $$

The color criterion becomes almost flat at $(g' - i')_0 = 0.779$ beyond $R_{maj} = 20'$. It is worth mentioning that the main radial gradient results are not dependent on the details of the color criterion and are little changed even if we adopt a single color at all radii, such as $(g - i)_0 = 0.8$ mag.

Figure 2(b) shows the $(g' - i')_0$ color of our GC samples as a function of major-axis distance from M87. The majority of the GCs are located within $R_{maj} = 65'$ where a strong sequence at $(g' - i')_0 \sim 0.7$ and a weak red tail are shown. Most of the red GCs are located within $R_{maj} = 40'$, while the blue GCs extend out to $R_{maj} = 180'$. We have 509 blue and 183 red GCs in our sample with the color criteria described above.

### 3. Age, [Z/H], [Fe/H], and [α/Fe] Measurements

The colors of GCs are an efficient tool to measure their metallicities. However, there are inherent uncertainties because GC colors also depend on GC ages. A spectral analysis has the advantage of measuring ages and metallicities of GCs more independently.

Lick line indices have been commonly used to measure ages and metallicities of old simple stellar populations with their integrated spectra. They are defined as the strengths of 25 prominent absorption lines seen in optical spectra of old stellar systems, established and refined by Burstein et al. (1984), Worthey et al. (1994), Worthey & Ottaviani (1997), and Trager et al. (1998).

We used the lick_ew code in the EZ_Ages package (Graves & Schiavon 2008) to measure the Lick indices from our spectra, adopting the passbands defined by Worthey & Ottaviani (1997) and Trager et al. (1998). This code smooths spectra to the resolution of the Lick system and estimates index errors based on error spectra (Cardiel et al. 1998).

We measured the mean ages, metallicities, and α-element abundances of the GC subgroups using a $\chi^2$ minimization technique (Proctor et al. 2004). This technique is based on the residuals between the Lick index measurements and the prediction from stellar population models. We used the flux-calibrated stellar population models of Lick indices presented by Thomas et al. (2011). The models were constructed with the ages ranging from 0.1 to 15 Gyr, the metallicities [Z/H] from -2.25 to +0.67, and the α-element abundances, [α/Fe], from -0.3 to +0.5. We measured the iron abundances by adopting the relation from Thomas et al. (2003), [Z/H] = [Fe/H] + 0.94 [α/Fe].

We used all Lick indices except $CN_1$, $CN_2$, Ca4227, and NaD to calculate the $\chi^2$ values. The stellar population model we adopted is calibrated with Galactic GCs that show an anomaly in the nitrogen abundance, which might not be applicable to typical simple stellar populations. The $CN_1$, $CN_2$, and Ca4227 are sensitive to nitrogen abundances, so they were excluded from the fitting. The NaD index is also excluded from
<figure><img src="image_6.png" /><figcaption>Figure 3. Examples of flux-calibrated spectra of (a) a blue GC and (d) a red GC that have S/N ~ 9. Middle and right panels show co-added spectra of GCs with S/Ns of ~30 and 80, respectively. The spectra were stacked for (b) three blue GCs at 2' < R<sub>maj</sub> < 2.4, (c) 53 blue GCs at 2' < R<sub>maj</sub> < 7.1, (e) two red GCs at 1.6 < R<sub>maj</sub> < 2.2, and (f) 45 red GCs at 1.6 < R<sub>maj</sub> < 9.2. All spectra were normalized at 5450-5550 Å.</figcaption></figure>

the fitting because it is strongly affected by interstellar absorption. We determined the $\chi^2$ values with the remaining Lick indices, using iterative $2\sigma$ clipping. Finally, 14–21 line indices of each co-added spectrum are used for the fitting.

Because our individual GC spectra typically do not have the required S/N for a robust stellar population analysis, we performed our analysis on co-added spectra. We stacked the spectra of various GC subgroups to obtain mean ages, metallicities, and $\alpha$-element abundances by taking the error-weighted mean values for each wavelength bin after strong sky emission line regions were masked. All spectra were de-redshifted prior to the co-addition. Figures 3(a) and (d) show flux-calibrated spectra of typical blue and red GC spectra, with S/N ~ 9 Å⁻¹ at 5000 Å, which is the median S/N for our sample. Panels (b), (c), (e), and (f) show stacked spectra of blue and red GCs with S/N of 30 and 80 as examples. The stacking analysis enables us to determine stellar population parameters from high-S/N spectra for GC subgroups. We have made the stacked spectra used in our analysis publicly available on Zenodo (10.5281/zenodo.6390374).

For the following analysis, the measurement errors for the ages, [Z/H], [Fe/H], and [α/Fe] were calculated using the bootstrap. For a given GC subgroup, we chose the same number of GCs in the group allowing replacement, stacked their spectra, and measured the stellar population parameters. We repeated this process 1000 times and took the 16th and 84th percentiles of the parameter distributions for lower and upper errors. The following results show the stellar population analysis for the GCs stacked using various criteria.

# 4. Results

## 4.1. Radial Gradients

Figures 4(a)–(d) show the mean age, [Z/H], [α/Fe], and [Fe/H] of the GC population in each radial bin as a function of major-axis distance from M87. The measurement values are listed in Table 1. The co-added spectrum of the GCs in each

radial bin has an S/N of 80. Note that the numbers of blue and red GCs in each radial bin are comparable within R = 14', but the blue GCs dominate in the outer region, where its number fraction rises up to ~80%.

The GC populations in all radial bins are older than 8 Gyr, but the ones in the outskirts have formal mean ages about 2 Gyr younger than those in the inner region. However, it is well known that changing horizontal branch morphology can affect age estimates of integrated stellar populations (e.g., Lee et al. 2000; Koleva et al. 2008; Cabrera-Ziri & Conroy 2022), especially for metal-poor GCs. Tests with our integrated spectra show that age estimates using only Hδ (and masking all other Balmer lines) are 1–2 Gyr younger than those using Hβ alone for the GC bins at large radii, which is a sign that blue horizontal branch stars may be having a moderate effect on age estimates (Schiavon et al. 2004). Therefore, although we report this apparent gradient in age, we urge some caution in their interpretation.

We find a steep radial gradient in the mean metallicities of the GCs within a certain radius and a flat gradient beyond. We derive broken linear fits for these data with the form as follows:

$$ [Z/H] = a + b \times \log\left(\frac{R_{\text{maj}}}{R_0}\right), \quad \text{if } R_{\text{maj}} < R_0, \\ = a, \quad \text{if } R_{\text{maj}} \ge R_0. \qquad (3) $$

The coefficient a is the value of [Z/H] at and beyond the scale distance R₀ (in kpc), and the coefficient b is the slope within R₀. The best-fit values we estimated are listed in Table 4. The entire GC sample shows a steep metallicity gradient with a slope of b = -0.84+0.15-0.11 within Rmaj ~ 22' (109 kpc) and a constant metallicity of [Z/H] = -1.34 beyond. This value of b means that, within R₀, [Z/H] changes by -0.84 dex for a factor of 10 in radius. The break point of Rmaj ~ 109 kpc is consistent with the distance at which the radial velocity dispersion of the GC system rises, found in Longobardi et al. (2018b), where the IGC population becomes dominant.
<figure><img src="image_7.png" /><figcaption>**Figure 4.** (a) Mean age, (b) mean [Z/H], (c) mean [α/Fe], and (d) mean [Fe/H] of the GC population as a function of major-axis distance from M87. Panels (e)–(h) show the same as panels (a)–(d), respectively, but for GC subpopulations: blue GCs (blue squares) and red GCs (red triangles). The solid lines in panels (b)–(d) and (f)–(h) represent the broken linear or linear fits for the GCs, blue and red GCs. The error bars on the y-axis represent the 1σ uncertainties from bootstrapping, while the errors bars in R<sub>maj</sub> along the x-axis represent the radial extent of each bin.</figcaption></figure>

The inner radial bins have estimated ages that saturate at the
limit of the models (15 Gyr). This model age zero-point issue
could potentially affect our metallicity measurements. To test
this, we artificially limit the model ages to be below 10 Gyr,
finding that the estimated [Z/H] increases by only ~0.1 dex.
This moderate increase is likely representative of the potential
effects of age saturation and is much smaller than the
magnitude of the detected gradient.

We also confirmed that the metallicities of individual GCs
with S/N > 20 are derived down to [Z/H] ~ -2 with
uncertainties of ~0.1. In addition, integrated metallicities of
the Milky Way GCs, measured using Schiavon et al. (2005)
spectra with the same method, are distributed below

[Z/H] = -1.6 and show a good agreement with the values
from the literature. We therefore conclude that the flattening of
the metallicity gradient in Figure 4(b) is real and not due to a
lack of metallicity sensitivity.

The α-element abundances also show a decreasing trend out to R<sub>maj</sub> ~ 22' and then become constant at larger distances. We fitted the data with the same functional form as used for [Z/H] (Equation (3)), with R<sub>0</sub> fixed to the value derived for [Z/H]. The coefficients *a* and *b* are derived to be *a* = 0.15 ± 0.04 and *b* = -0.34 ± 0.11. Described in words, this value of *b* reflects the fact that GCs in M87's inner region are roughly twice as abundant in alpha-elements relative to iron compared to GCs at and beyond R<sub>0</sub>. We perform a similar analysis for the radial
<table><caption>**Table 1**
Mean Ages, [Z/H], [α/Fe], and [Fe/H] for the Entire GC Sample</caption><thead><tr><th>R<sub>maj</sub> Range (kpc)</th><th>R&#773;<sub>maj</sub><sup>a</sup> (kpc)</th><th>Age (Gyr)</th><th>[Z/H] (dex)</th><th>[α/Fe] (dex)</th><th>[Fe/H] (dex)</th></tr></thead><tbody><tr><td>8–18</td><td>14</td><td>15.0<sup>+0.1</sup><sub>-0.1</sub></td><td>-0.55<sup>+0.11</sup><sub>-0.12</sub></td><td>0.41<sup>+0.09</sup><sub>-0.14</sub></td><td>-0.94<sup>+0.10</sup><sub>-0.05</sub></td></tr><tr><td>17–38</td><td>28</td><td>15.0<sup>+0.1</sup><sub>-0.1</sub></td><td>-0.95<sup>+0.08</sup><sub>-0.03</sub></td><td>0.42<sup>+0.07</sup><sub>-0.05</sub></td><td>-1.34<sup>+0.09</sup><sub>-0.05</sub></td></tr><tr><td>38–58</td><td>48</td><td>15.0<sup>+0.1</sup><sub>-1.1</sub></td><td>-1.04<sup>+0.08</sup><sub>-0.05</sub></td><td>0.29<sup>+0.05</sup><sub>-0.08</sub></td><td>-1.31<sup>+0.12</sup><sub>-0.07</sub></td></tr><tr><td>58–75</td><td>65</td><td>10.2<sup>+4.8</sup><sub>-1.0</sub></td><td>-1.01<sup>+0.07</sup><sub>-0.14</sub></td><td>0.09<sup>+0.18</sup><sub>-0.04</sub></td><td>-1.09<sup>+0.07</sup><sub>-0.28</sub></td></tr><tr><td>75–90</td><td>83</td><td>12.6<sup>+2.4</sup><sub>-1.3</sub></td><td>-1.20<sup>+0.07</sup><sub>-0.10</sub></td><td>0.27<sup>+0.11</sup><sub>-0.06</sub></td><td>-1.45<sup>+0.08</sup><sub>-0.15</sub></td></tr><tr><td>90–108</td><td>99</td><td>13.0<sup>+1.0</sup><sub>-1.9</sub></td><td>-1.40<sup>+0.11</sup><sub>-0.01</sub></td><td>0.27<sup>+0.09</sup><sub>-0.15</sub></td><td>-1.65<sup>+0.20</sup><sub>-0.04</sub></td></tr><tr><td>109–139</td><td>126</td><td>11.6<sup>+3.4</sup><sub>-1.7</sub></td><td>-1.29<sup>+0.08</sup><sub>-0.15</sub></td><td>0.07<sup>+0.10</sup><sub>-0.13</sub></td><td>-1.36<sup>+0.12</sup><sub>-0.17</sub></td></tr><tr><td>139–175</td><td>154</td><td>11.4<sup>+3.6</sup><sub>-1.1</sub></td><td>-1.32<sup>+0.05</sup><sub>-0.13</sub></td><td>0.14<sup>+0.10</sup><sub>-0.11</sub></td><td>-1.45<sup>+0.08</sup><sub>-0.14</sub></td></tr><tr><td>175–232</td><td>197</td><td>9.8<sup>+5.2</sup><sub>-0.4</sub></td><td>-1.33<sup>+0.04</sup><sub>-0.23</sub></td><td>0.07<sup>+0.10</sup><sub>-0.12</sub></td><td>-1.40<sup>+0.04</sup><sub>-0.22</sub></td></tr><tr><td>234–376</td><td>279</td><td>11.7<sup>+2.3</sup><sub>-0.8</sub></td><td>-1.35<sup>+0.01</sup><sub>-0.12</sub></td><td>0.17<sup>+0.11</sup><sub>-0.15</sub></td><td>-1.51<sup>+0.10</sup><sub>-0.14</sub></td></tr><tr><td>381–837</td><td>598</td><td>12.5<sup>+1.5</sup><sub>-0.9</sub></td><td>-1.40<sup>+0.06</sup><sub>-0.09</sub></td><td>0.18<sup>+0.10</sup><sub>-0.13</sub></td><td>-1.57<sup>+0.13</sup><sub>-0.11</sub></td></tr></tbody></table>

**Note.**

<table><caption>**Table 2**
Mean Ages, [Z/H], [α/Fe], and [Fe/H] for the Blue GC Sample</caption><thead><tr><th>R<sub>maj</sub> Range (kpc)</th><th>R&#773;<sub>maj</sub><sup>a</sup> (kpc)</th><th>Age (Gyr)</th><th>[Z/H] (dex)</th><th>[α/Fe] (dex)</th><th>[Fe/H] (dex)</th></tr></thead><tbody><tr><td>9–25</td><td>19</td><td>15.0<sup>+0.1</sup><sub>-0.7</sub></td><td>-0.95<sup>+0.12</sup><sub>-0.04</sub></td><td>0.38<sup>+0.12</sup><sub>-0.03</sub></td><td>-1.31<sup>+0.08</sup><sub>-0.03</sub></td></tr><tr><td>25–54</td><td>40</td><td>14.0<sup>+1.0</sup><sub>-0.9</sub></td><td>-1.25<sup>+0.07</sup><sub>-0.03</sub></td><td>0.31<sup>+0.03</sup><sub>-0.14</sub></td><td>-1.54<sup>+0.18</sup><sub>-0.14</sub></td></tr><tr><td>54–77</td><td>65</td><td>10.3<sup>+2.7</sup><sub>-0.1</sub></td><td>-1.28<sup>+0.01</sup><sub>-0.18</sub></td><td>0.14<sup>+0.12</sup><sub>-0.16</sub></td><td>-1.41<sup>+0.08</sup><sub>-0.16</sub></td></tr><tr><td>77–94</td><td>85</td><td>13.5<sup>+0.5</sup><sub>-1.7</sub></td><td>-1.44<sup>+0.09</sup><sub>-0.04</sub></td><td>0.18<sup>+0.12</sup><sub>-0.16</sub></td><td>-1.61<sup>+0.19</sup><sub>-0.16</sub></td></tr><tr><td>94–112</td><td>103</td><td>9.9<sup>+2.6</sup><sub>-0.5</sub></td><td>-1.35<sup>+0.01</sup><sub>-0.18</sub></td><td>0.28<sup>+0.13</sup><sub>-0.15</sub></td><td>-1.61<sup>+0.09</sup><sub>-0.15</sub></td></tr><tr><td>112–142</td><td>128</td><td>13.0<sup>+1.0</sup><sub>-1.0</sub></td><td>-1.48<sup>+0.07</sup><sub>-0.05</sub></td><td>-0.10<sup>+0.14</sup><sub>-0.15</sub></td><td>-1.39<sup>+0.14</sup><sub>-0.15</sub></td></tr><tr><td>141–181</td><td>158</td><td>12.9<sup>+0.1</sup><sub>-1.5</sub></td><td>-1.57<sup>+0.10</sup><sub>-0.03</sub></td><td>0.05<sup>+0.19</sup><sub>-0.18</sub></td><td>-1.62<sup>+0.18</sup><sub>-0.18</sub></td></tr><tr><td>181–241</td><td>206</td><td>10.8<sup>+3.2</sup><sub>-0.9</sub></td><td>-1.47<sup>+0.05</sup><sub>-0.15</sub></td><td>0.19<sup>+0.06</sup><sub>-0.24</sub></td><td>-1.65<sup>+0.19</sup><sub>-0.24</sub></td></tr><tr><td>241–421</td><td>301</td><td>15.0<sup>+0.1</sup><sub>-5.0</sub></td><td>-1.61<sup>+0.15</sup><sub>-0.05</sub></td><td>0.15<sup>+0.08</sup><sub>-0.30</sub></td><td>-1.75<sup>+0.33</sup><sub>-0.30</sub></td></tr><tr><td>423–829</td><td>615</td><td>12.6<sup>+0.4</sup><sub>-0.9</sub></td><td>-1.59<sup>+0.09</sup><sub>-0.06</sub></td><td>0.09<sup>+0.20</sup><sub>-0.13</sub></td><td>-1.67<sup>+0.14</sup><sub>-0.13</sub></td></tr></tbody></table>

**Note.**

<table><caption>**Table 3**
Mean Ages, [Z/H], [α/Fe], and [Fe/H] for the Red GC Sample</caption><thead><tr><th>R<sub>maj</sub> Range (kpc)</th><th>R&#773;<sub>maj</sub><sup>a</sup> (kpc)</th><th>Age (Gyr)</th><th>[Z/H] (dex)</th><th>[α/Fe] (dex)</th><th>[Fe/H] (dex)</th></tr></thead><tbody><tr><td>8–29</td><td>18</td><td>15.0<sup>+0.1</sup><sub>-0.1</sub></td><td>-0.52<sup>+0.27</sup><sub>-0.01</sub></td><td>0.24<sup>+0.26</sup><sub>-0.01</sub></td><td>-0.75<sup>+0.05</sup><sub>-0.06</sub></td></tr><tr><td>30–62</td><td>48</td><td>15.0<sup>+0.1</sup><sub>-1.0</sub></td><td>-0.72<sup>+0.11</sup><sub>-0.05</sub></td><td>0.33<sup>+0.07</sup><sub>-0.06</sub></td><td>-1.03<sup>+0.11</sup><sub>-0.06</sub></td></tr><tr><td>61–103</td><td>79</td><td>15.0<sup>+0.1</sup><sub>-3.0</sub></td><td>-0.72<sup>+0.10</sup><sub>-0.04</sub></td><td>0.34<sup>+0.04</sup><sub>-0.06</sub></td><td>-1.04<sup>+0.12</sup><sub>-0.04</sub></td></tr><tr><td>104–253</td><td>165</td><td>9.2<sup>+1.7</sup><sub>-2.2</sub></td><td>-0.67<sup>+0.12</sup><sub>-0.11</sub></td><td>0.26<sup>+0.04</sup><sub>-0.07</sub></td><td>-0.91<sup>+0.15</sup><sub>-0.12</sub></td></tr><tr><td>259–837</td><td>488</td><td>15.0<sup>+0.1</sup><sub>-3.7</sub></td><td>-0.84<sup>+0.17</sup><sub>-0.04</sub></td><td>0.25<sup>+0.06</sup><sub>-0.11</sub></td><td>-1.07<sup>+0.23</sup><sub>-0.05</sub></td></tr></tbody></table>

**Note.**

trend of [Fe/H]. Here we also find a significant gradient within $R_0$. Although we again fixed $R_0$ to 109 kpc, allowing $R_0$ to vary freely produces a very similar value. The coefficients *a* and *b* for [Fe/H] are estimated to be $a = -1.48^{+0.03}_{-0.07}$ and $b = -0.52^{+0.06}_{-0.14}$. The magnitudes of the slopes for [α/Fe] and [Fe/H] are comparable within the uncertainties, showing that they each contribute roughly equally to the [Z/H] gradient.

We also investigate the radial trends of the stellar populations of GC subpopulations where the co-added spectrum in each radial bin has an S/N of 70. The mean ages, [Z/H], [α/Fe], and [Fe/H] for blue and red GC samples are listed in Tables 2 and 3, respectively. The blue and red GCs in most radial bins have mean ages older than 10 Gyr

(Figure 4(e)). Figure 4(f) shows the metallicity gradients of the blue and red GCs. The red GCs are clearly more metal-rich than the blue GCs, which indicates that the $(g' - i')$ color corresponds well to the average metallicity.

We fit linear relations for the metallicity and α-element abundance gradients of the blue subpopulation in the same way as they were derived for the entire sample. For the red GCs, we do a single linear fit, as there is no convincing evidence of a break point. The best-fit coefficients are listed in Table 4. The blue GCs have a steep total metallicity gradient with a slope of $b = -0.61^{+0.05}_{-0.21}$ within $R_{maj} \sim 34'$ (165 kpc) and a constant metallicity of $[Z/H] = -1.56$ beyond, which is a similar trend to that of the entire sample. On the other hand, the red GCs

<footnote>a Luminosity-weighted mean major-axis distance of GCs in each radial bin.</footnote>

<footnote>a Luminosity-weighted mean major-axis distance of blue GCs in each radial bin.</footnote>

<footnote>a Luminosity-weighted mean major-axis distance of red GCs in each radial bin.</footnote>
<table><caption>**Table 4**
Linear Fit Results for Radial Gradients of Mean [Z/H], [α/Fe], and [Fe/H] for the GC Subsample</caption><thead><tr><th>Sample</th><th>R<sub>maj</sub> Range (kpc)</th><th>R<sub>0</sub> (kpc)</th><th>a<sub>[Z/H]</sub></th><th>b<sub>[Z/H]</sub></th><th>a<sub>[α/Fe]</sub></th><th>b<sub>[α/Fe]</sub></th><th>a<sub>[Fe/H]</sub></th><th>b<sub>[Fe/H]</sub></th></tr></thead><tbody><tr><td rowspan="2">All</td><td>8 – R<sub>0</sub></td><td>109<sup>+82</sup><sub>-2</sub></td><td>-1.34<sup>+<0.01</sup><sub>-0.12</sub></td><td>-0.84<sup>+0.15</sup><sub>-0.11</sub></td><td>0.15<sup>+0.04</sup><sub>-0.04</sub></td><td>-0.34<sup>+0.11</sup><sub>-0.11</sub></td><td>-1.48<sup>+0.03</sup><sub>-0.07</sub></td><td>-0.52<sup>+0.06</sup><sub>-0.14</sub></td></tr><tr><td>R<sub>0</sub> – 837</td><td>109<sup>+82</sup><sub>-2</sub></td><td>-1.34<sup>+<0.01</sup><sub>-0.12</sub></td><td>0</td><td>0.15<sup>+0.04</sup><sub>-0.04</sub></td><td>0</td><td>-1.48<sup>+0.03</sup><sub>-0.07</sub></td><td>0</td></tr><tr><td rowspan="2">Blue</td><td>9 – R<sub>0</sub></td><td>165<sup>+32</sup><sub>-60</sub></td><td>-1.56<sup>+0.05</sup><sub>-0.05</sub></td><td>-0.61<sup>+0.05</sup><sub>-0.21</sub></td><td>0.09<sup>+0.04</sup><sub>-0.10</sub></td><td>-0.30<sup>+0.08</sup><sub>-0.17</sub></td><td>-1.64<sup>+0.09</sup><sub>-0.05</sub></td><td>-0.32<sup>+0.16</sup><sub>-0.10</sub></td></tr><tr><td>R<sub>0</sub> – 829</td><td>165<sup>+32</sup><sub>-60</sub></td><td>-1.56<sup>+0.05</sup><sub>-0.05</sub></td><td>0</td><td>0.09<sup>+0.04</sup><sub>-0.10</sub></td><td>0</td><td>-1.64<sup>+0.09</sup><sub>-0.05</sub></td><td>0</td></tr><tr><td>Red</td><td>8–837</td><td>100<sup>a</sup></td><td>-0.70<sup>+0.09</sup><sub>-0.01</sub></td><td>-0.18<sup>+0.09</sup><sub>-0.15</sub></td><td>0.28<sup>+0.04</sup><sub>-0.04</sub></td><td>-0.01<sup>+0.05</sup><sub>-0.18</sub></td><td>-0.97<sup>+0.08</sup><sub>-0.02</sub></td><td>-0.16<sup>+0.15</sup><sub>-0.05</sub></td></tr></tbody></table>

Note. The functional form of linear fits is $a + b \times \log(R_{maj}/R_0)$.

a For red GCs, R₀ is not a break point of segmented linear fits, but simply a reference point.

have a much shallower gradient in [Z/H], with a slope of $b = -0.18^{+0.09}_{-0.15}$. A value of [Z/H] ≈ -0.7 describes the red GCs across a large range in radius.

Figure 4(g) shows the radial trends of α-element-to-iron abundance ratios of the blue and red GCs. The blue GCs show a noticeable gradient with a slope of $b = -0.30^{+0.08}_{-0.17}$ within R<sub>maj</sub> ~ 34' (165 kpc), while the red GCs have a constant [α/Fe] value of ~0.3 in all radial ranges with a radial gradient slope of $b = -0.01^{+0.05}_{-0.18}$. The mean [α/Fe] value of the outermost blue GCs is +0.09, similar to the solar abundance.

Recently, Villaume et al. (2020) presented a stellar population study of M87 GCs based on full spectrum fitting of individual GC spectra. They found no significant radial gradient in [Fe/H] for either blue or red GCs. A direct comparison to this work is complicated by the fact that the sample of Villaume et al. (2020) covers a smaller range in radius ($R < 150$ kpc) and a gap in their data at 50 kpc that motivated their choice of break point in their fitting. Our [Fe/H] gradient for the blue GCs (-0.32 dex per decade in radius) is consistent with their findings, within the uncertainties, although a more homogenous analysis between the two data sets is certainly desirable in the future. Villaume et al. (2020) also find a strong negative gradient in the blue GCs for certain α-elements, like [Mg/Fe], which is similar to what we see in this work.

## 4.2. Azimuthal Variations

The Virgo Cluster core around M87 is known to be a place of active galaxy and cluster assembly that is likely not well mixed (e.g., Mihos et al. 2017). Our wide spatial coverage of GCs in this region allows us to investigate the possibility that there may also be azimuthal variations of mean ages and total metallicities of the GC population in addition to radial trends. We take ellipses centered on M87 that have a position angle of 155° and an ellipticity of 0.4 (Ferrarese et al. 2006; Janowiecki et al. 2010). We divide these ellipses into four sectors azimuthally: northeast (NE), northwest (NW), southeast (SE), and southwest (SW) (Figure 5). We further subdivide these sectors into four regions along the major-axis distance ($R_{maj} < 10'$, $10' < R_{maj} < 25'$, $25' < R_{maj} < 65'$, and $R_{maj} > 65'$). Figure 5 shows the S/N of the co-added spectra of the blue and red GCs in subregions, ranging from 6 to 68. We measured the mean ages and total metallicities for 16 and 14 co-added spectra for the blue and red GCs, respectively.

The blue and red GCs in 26 out of 30 subregions are older than 10 Gyr, not showing specific trends. In addition to the sectors, we investigated the substructure that Powalka et al. (2018) identified as the region that relatively young GCs

<figure><img src="image_8.png" /><figcaption>**Figure 5.** (a) S/N of the co-added spectra for the blue GC population in spatial bins. The solid-line ellipses are centered on M87 and have semimajor axes of 10', 25', and 65'; a position angle of 155°; and an ellipticity of 0.4. The dashed-line and dotted-line ellipses are centered on M86 and M84 with semimajor axes of 23' and 15', corresponding to 10 effective radii of M86 and M84, respectively. The position angle and ellipticity of M86 are 123° and 0.21, and those of M84 are 128° and 0.14 (Ferrarese et al. 2020). The radial solid lines indicate the position angle criteria used to divide the survey region into different sectors. The numbers in each subregion indicate the numbers of the co-added GCs. (b) Same as panel (a), but for the red GC population.</figcaption></figure>

populate, located at 11'5 south of M87 with a diameter of ~6'. We found 41 GCs in that region and measured their mean age with the co-added spectrum. They have a mean age of 9 Gyr, which is indeed a few Gyr smaller than the azimuthally averaged age of other GCs at the same distance from M87. Further kinematic studies of the GCs in this region will be helpful to reveal the origin of any younger GC population.

The two GC populations show radial gradients in metallicity for each sector. Figure 6 shows the radial trends in metallicity of the blue and red GCs in four sectors. The radial trends of the
<figure><img src="image_9.png" /></figure>

**Figure 6.** [Z/H] vs. major-axis distance from M87 for the blue GCs (squares) and red GCs (triangles) in four different sectors: (a) NE, (b) NW, (c) SE, and (d) SW. The blue and red dashed lines are the same in all plots and indicate the radial gradients in the metallicities of the full blue and red GC subpopulations, respectively (Figure 4(f)). The filled and open symbols represent the GC samples with co-added spectral S/N higher and lower than 30, respectively.

blue GCs in all sectors are consistent with those derived with the entire blue GC sample (Figure 4(e)). The red GCs also show shallow radial gradients in general, but the red GCs in the NW region of M87—the direction of M86 and M84—have slightly elevated mean metallicities compared to those in other sectors (Figure 6(b)). There is a possibility that this metal-rich population originates from massive galaxies in the northwest region of M87. However, the difference is not significant enough for us to draw a definite conclusion.

# 5. Discussion: The Assembly History of M87 and the Virgo Cluster Core

We have estimated ages and metallicities from co-added spectra of GC subpopulations and have found several notable features, especially in GC metallicity variations as a function of projected distance from M87. The GCs in M87 and the intracluster region are mostly old, potentially allowing us to trace the mass assembly history of M87 to early epochs (>10 Gyr). In this section, we discuss the influence of two different parameters related to the assembly history of the Virgo core region and M87 itself: the masses of GC progenitor host galaxies and the dependence on environment.

## 5.1. Gradients and GC Progenitor Galaxy Masses

The chemical abundances of stars and star clusters are intimately tied to the mass of the galaxy within which they are formed. The stellar mass–metallicity relation of galaxies has connected these quantities since its discovery (Lequeux et al. 1979), and the relation was confirmed with gas-phase metallicities, especially for star-forming galaxies (e.g., Tremonti et al. 2004; Mannucci et al. 2010; Andrews & Martini 2013). A series of studies have also clarified and extended the relation

with stellar metallicities for Sloan Digital Sky Survey star-forming galaxies (Gallazzi et al. 2005) down to Local Group dwarf spheroidal and irregular galaxies (Kirby et al. 2013). Although the metallicities of the GCs are about 0.8–1 dex lower than those of the stars in a galaxy (Jordán et al. 2004; Lamers et al. 2017), we expect that the metallicities of the GCs also represent masses of progenitors where the GCs were formed. Therefore, one interpretation of the metallicity radial gradients seen in the blue and red GCs is that both blue and red GCs in the inner region of M87 were formed in more massive galaxies than those in the outskirts.

If the GCs in massive galaxies like M87 largely originate from merged and accreted galaxies that were disrupted (e.g., Côté et al. 1998) and GC metallicity correlates with progenitor host mass, then a metallicity gradient suggests a radial gradient in the stellar masses of disrupted progenitors. This is a natural consequence of the efficiency of dynamical friction as a function of the merger ratio (Amorisco 2019). Merged galaxies will spiral into the center owing to dynamical friction until they are disrupted. Not only do more massive progenitors experience dynamical friction more efficiently, but they can also resist complete disruption until they are deeper in the potential well of the larger galaxy. In this way, one expects more metal-rich stars and star clusters to be deposited closer to the center of the accreting galaxy, creating a metallicity gradient.

To check whether this explanation can work in practice, we use the Illustris cosmological simulation (Vogelsberger et al. 2014a, 2014b), where we analyze the GC system of a simulated galaxy cluster with a virial mass consistent with that of the Virgo Cluster as an example. Ramos-Almendares et al. (2020) tagged GCs to individual galaxies in the simulated galaxy at their infall time, following a Hernquist density profile (Hernquist 1990). The infall time here is defined as the epoch when a given galaxy can no longer be identified as the central galaxy of its own friends-of-friends group. They adopted different values of the characteristic density and scale length in the Hernquist profile for blue and red GC systems. The relation between the total mass in GCs and host galaxy mass (Harris et al. 2015; Hudson et al. 2015) is taken into account to assign masses to GCs. The detailed GC tagging analysis is described in Ramos-Almendares et al. (2020). In our analysis, we only study “orphan” GCs—GCs for which the progenitor host galaxies are destroyed at the present day.

We investigate where orphan GCs are located as a function of their progenitors’ masses. The mean magnitude of our spectroscopic sample is $\langle g' \rangle = 21.5$ mag, with 99% of the sample having $g' < 22.5$ mag. The latter limit corresponds to a mass of $M = 4 \times 10^5 M_{\odot}$ with an assumption of a GC mass-to-light ratio in the g-band of 2. We therefore only select the simulated galaxies with GC system masses higher than $4 \times 10^5 M_{\odot}$ for a fair comparison. These galaxies have stellar masses ranging from $3 \times 10^8 M_{\odot}$ to $10^{11} M_{\odot}$. Assuming that the mass of a GC is $4 \times 10^5 M_{\odot}$, we randomly chose GC particles until their total mass was equal to the total mass of GCs in each galaxy. We repeated this particle sampling 100 times to determine the mean trend of progenitors’ masses as a function of projected distance from the central galaxy.

Figure 7 shows the mean stellar masses at infall time of the progenitor galaxies of orphan GCs as a function of their present-day projected clustercentric distance. We find a steep gradient in the mean progenitor masses for blue orphan GCs
<figure><img src="image_10.png" /><figcaption>**Figure 7.** (a) Mean stellar masses of progenitors at infall time as a function of projected clustercentric distance for blue orphan GCs in a simulated galaxy cluster (Ramos-Almendares et al. 2020). (b) Same as panel (a), but for red orphan GCs. The solid lines indicate the running averages of the progenitor masses with moving bins of *N* = 100. The shaded regions represent the standard deviation from the running averages. Note that this plot does not include any GC populations formed in situ.</figcaption></figure>

located at 50 kpc $< R_{\text{proj}} < 300$ kpc, ranging from $10^{10.5}$ to $10^{9.3} M_{\odot}$. This implies that the blue orphan GCs in the inner region of the central galaxy originate from more massive satellites that are more likely to sink further into the center than those in the outer region. We cannot constrain the gradient in the innermost region ($R_{\text{proj}} < 50$ kpc) because in this analysis we do not include the “in situ” population that has formed at the center and would be the most metal-rich population in the central galaxy. On the other hand, the red orphan GCs do not show a steep gradient in the mean progenitor masses because the progenitor masses of the red GCs have a much narrower range than those of the blue GCs. The red GCs are rarely formed in galaxies with $M_{* \text{, infall}} < 10^{10} M_{\odot}$. This trend matches well with the observed radial gradient in metallicity for the red GCs being shallower than that of the blue GCs.

In addition, we also found that the orphan GCs, especially blue ones, from progenitors with the earliest infall time ($z \sim 3.6$) are located at the innermost region in this simulation, while those from recent infallers with a mean infall time of $z \sim 2.5$ are distributed in the outer region. This hints at the observed age difference of a few Gyr between the inner and outer GCs.

While the simulations have their limitations, this example illustrates how including GCs in cosmological galaxy simulations, in a way that allows for studies of GC properties as a function of their spatial distributions (e.g., Ramos-Almendares et al. 2020; Reina-Campos et al. 2021; Chen & Gnedin 2022), makes available new diagnostics for galaxy assembly history.

## 5.2. GC Occupation Fraction and Gradient Flattening

In addition to the steeper radial metallicity gradient for blue GCs, we have found that the gradient flattens beyond $\sim 165$ kpc (Figure 4(f)). Since the mass–metallicity relation extends to quite low metallicities and masses, and the gravitational

processes that govern dynamical friction and galaxy disruption are also continuous, it might seem surprising that we should find any break in these radial trends. We show below that this could be related to the GC occupation fraction—the fraction of galaxies that contain GCs—as a function of galaxy mass.

Sánchez-Janssen et al. (2019) presented the star cluster occupation fraction of Virgo early-type galaxies as a function of galaxy stellar mass (see their Figure 6). All galaxies with $M_{*} > 10^{9} M_{\odot}$ host GCs with $g' < 24.5$ mag (roughly the mean of the GC luminosity function), while only half of the galaxies with $M_{*} \sim 10^{7} M_{\odot}$ do. Because of the flux limit of our survey, we only observed relatively bright GCs, with a mean sample magnitude of $\langle g' \rangle = 21.5$ mag, and with 99% of the sample having $g' < 22.5$ mag. The relevant GC occupation fraction for determining whether an accreted progenitor contributes GCs to our observed sample is therefore whether it contains a GC with $g' < 22.5$ mag. Using the GC data from the ACS Virgo Cluster Survey (Côté et al. 2004; Jordán et al. 2009), we determine the GC occupation fraction for Virgo galaxies as a function of GC magnitude and galaxy stellar mass. For galaxies with $M_{*} \sim 10^{10} M_{\odot}$, the occupation fraction for GCs with $g' < 22.5$ mag is $\sim 80\%$, but for galaxies with $M_{*} \sim 10^{9} M_{\odot}$, the occupation fraction drops to $\sim 25\%$. Galaxies with even lower stellar masses are then not expected to host many GCs with $g' < 22.5$ mag. Therefore, we expect that beyond $\sim 100–150$ kpc, while many lower-mass galaxies are being accreted and disrupted, only a small fraction of them can contribute GCs. The flattening of the gradients at large distances is a natural by-product of the fact that, below a certain mass, the low-mass galaxies that would host even more metal-poor GCs have few to no GCs. We see a sign of this in the simulated GC system shown in Figure 7(a). Beyond a certain distance, the mean progenitor masses of blue orphan GCs become constant because in this simulated galaxy cluster
the galaxies disrupted in the outskirts have lower masses and do not host any GCs more massive than $4 \times 10^5 M_\odot$. The location of the break point likely depends on the mass and assembly history of individual galaxy clusters. While the flux limit of our survey may affect the exact [Z/H] or [α/Fe] value at which the gradient flattens, we expect that the general trend would be similar even if one could observe fainter GCs, and this could be tested with future observations.

We present this as an alternative to a more traditional explanation for stellar population gradients, which is that they occur during dissipational star formation, where more rapid star formation and gas recycling toward the galaxy center cause more rapid chemical evolution. Evidence for steep metallicity gradients in “relic” galaxies with little apparent accreted populations (Martín-Navarro et al. 2015; Beasley et al. 2018) suggests that in situ star formation likely plays a role in these gradients, especially for more metal-rich populations. However, given the complexity of modeling a dissipational collapse, we find it interesting to show that the trends we see in the M87 GC system can be explained using a picture dominated by dry merging. This is especially relevant because it is expected that the intracluster light is dominated by stars acquired through dissipationless accretion (Purcell et al. 2007). It is likely, however, that dissipational collapse plays some role in forming the gradients we observe, especially in the inner regions, and it will be interesting to explore that possibility with future simulations.

### 5.3. [α/Fe] and Dependence on Environment

The [α/Fe] ratios allow us to understand the star formation timescales in the host environments. There is a well-known relation between α-element abundances and central velocity dispersion values for galaxies. The larger velocity dispersion the galaxies have, the higher [α/Fe] values they have (Thomas et al. 2005; Annibali et al. 2007; McDermid et al. 2015). However, for low-mass galaxies, the local density could be a more critical parameter than the mass. Liu et al. (2016) measured the α-element abundances of 11 low-mass dwarf galaxies in the Virgo Cluster and found that the [α/Fe] values of dwarf galaxies closer to M87 are higher than those distant from M87.

Similarly, we find moderate negative radial gradients in mean [α/Fe] values, especially for the blue GC sample, a population that likely originates from low-mass galaxies. This suggests that the blue GCs in the central region of M87 were formed in progenitors where rapid star formation occurred. Figure 8 shows the comparison of [α/Fe] gradients of GCs from this study and dwarf galaxies from Liu et al. (2016). The gradients are similar to each other, although the GCs are more centrally concentrated (in projected radius) than the Virgo dwarfs.

Following the simulation result above, the blue GCs in the inner region of M87 might originate from relatively massive galaxies that became satellites of the central galaxy at high redshift. Intense star formation and rapid quenching occurred in these massive progenitors in the high-z universe, resulting in high [α/Fe] abundances in their stars and star clusters. Similarly, dwarf galaxies in dense environments have enhanced α-element abundances (Chilingarian et al. 2008; Smith et al. 2009; Liu et al. 2016). There is a possibility that the most metal-rich of the blue GCs in the inner region originate from

<figure><img src="image_11.png" /><figcaption>**Figure 8.** [α/Fe] values of blue GCs (blue squares) and Virgo dwarf galaxies (green circles; Liu et al. 2016) as a function of major-axis distance from M87.</figcaption></figure>

dwarf galaxies that were quenched quickly by the dense environment (Mistani et al. 2016).

### 5.4. Local Substructures

The Virgo Cluster is a dynamically young galaxy cluster with several subclusters (Binggeli et al. 1985, 1987, 1993). Mihos et al. (2005, 2017) presented the deep imaging survey showing the diffuse intracluster light in the Virgo core region, indicating the ongoing interaction between galaxies. Related to this dynamical status, it is expected that the GCs can be stripped not only from low-mass dwarf galaxies but also from massive galaxies. The observed metal-rich GCs imply the ongoing interaction between M87 and nearby massive galaxies.

We detected a hint that the metal-rich GCs located in the NW region of M87 (Figure 6) are slightly more metal-rich than in other areas of the Virgo core. This direction is toward massive early-type galaxies, M86 and M84. We suspect that these metal-rich GCs might be associated with M84 because we exclude the GCs that have lower radial velocities corresponding to the systemic velocity of M86 (Figure 1(a)). We defer to a future paper an investigation of the many GCs associated with known substructures, particularly the M86 group.

## 6. Summary

We present a stellar population analysis of the GCs in the central region of the Virgo Cluster, especially around M87, using the MMT/Hectospec spectra. We select 692 GCs that belong to the main body of the Virgo Cluster to investigate the stellar populations of GCs, excluding the GCs associated with substructures. The GCs are divided into several groups with various criteria for colors, major-axis distances, and position angles. We measured mean ages, total metallicities, [Fe/H], and α-element abundances of the GC subpopulations using their co-added spectra based on Lick indices. The main findings of this study are summarized as follows.
1. The GCs as a whole follow a two-component radial metallicity ([Z/H]) trend, with a steep gradient following a power-law slope of $-0.84_{-0.11}^{+0.15}$ within $R_{maj} = 22'$ (109 kpc), and a flat gradient with a mean metallicity of $[Z/H] \sim -1.34$ beyond. The trends for [Fe/H] and $[\alpha/Fe]$ are similar, with inner region slopes of $-0.52_{-0.14}^{+0.06}$ and $-0.34 \pm 0.11$ and outer flattenings at [Fe/H] = $-1.48_{-0.07}^{+0.03}$ and $[\alpha/Fe] = 0.15 \pm 0.04$, respectively.

2. When analyzing gradients in GC color subpopulations (blue and red), we find that the two-component trend is clearly seen only in the blue GC system, while the red GCs are best fit by a single component. For [Z/H], the blue GCs exhibit an inner region power-law slope of $-0.61_{-0.21}^{+0.05}$ within $R_{maj} = 34'$ (165 kpc) and a mean metallicity of [Z/H] = $-1.56 \pm 0.05$ beyond. The red GCs show a much shallower radial gradient in [Z/H] than do the blue GCs, with a power-law slope of $-0.18_{-0.15}^{+0.09}$ and no obvious flattening at large radius.

3. The [Fe/H] and $[\alpha/Fe]$ gradients also behave differently for the blue and red GC subpopulations. The blue GCs follow the same two-component trend, with inner region slopes of $-0.32_{-0.10}^{+0.16}$ and $-0.30_{-0.17}^{+0.08}$ and outer flattenings at [Fe/H] = $-1.64_{-0.05}^{+0.09}$ and $[\alpha/Fe] = +0.09_{-0.10}^{+0.04}$, respectively. These results show that the [Z/H] gradient seen in the blue GCs is driven in roughly equal portions by [Fe/H] and $[\alpha/Fe]$. By contrast, the red GCs exhibit only a weak gradient in [Fe/H] (slope of $-0.16_{-0.05}^{+0.15}$) and no gradient in $[\alpha/Fe]$.

4. We investigated the azimuthal variation of GC metalli-
cities and detected a hint of local substructures toward the
massive early-type galaxy M84, consisting of metal-
rich GCs.

5. The GCs in the Virgo core are mostly old (>8 Gyr) regardless of their colors and locations. The outer GCs are formally about 2 Gyr younger than inner GCs, but there is a possibility that the presence of blue horizontal branch stars is causing this apparent age gradient.

We interpreted the radial gradient of GCs in metallicity and
$[\alpha/Fe]$ in the context of the assembly history of M87 and the
Virgo Cluster core. We used the Illustris cosmological
simulations to investigate masses of disrupted GC progenitors
in a simulated cluster.

1. According to the mass-metallicity relation of galaxies, the more metal-rich blue GCs at smaller galactocentric distances were probably formed in more massive galaxies. We expect that these more massive progenitors were better able to survive disruption as they approached the galaxy center due to dynamical friction. From the simulation, we found that the GCs from massive progenitors sink further into the central galaxy, while GCs from lower-mass progenitors are located preferentially at larger distances.

2. The flattening of metallicity gradients of blue GCs at large distances ($R \gtrsim 100$ kpc) is a natural consequence if the lower-mass dwarfs that could contribute GCs with even lower metallicity actually host few or no GCs.

3. We suggested that the dense environment around M87 could result in the high $[\alpha/Fe]$ values of the GCs in the central region, mirroring the environmental dependence of the $[\alpha/Fe]$ of Virgo dwarf galaxies (Liu et al. 2016).

Overall, we suggest a scenario where more massive progenitors accreted into the central galaxy (M87) contribute a significant fraction of the GCs in the inner region. Future simulation studies on various galaxy clusters will be helpful to understand the relation between the assembly history of a host galaxy cluster and chemical properties of the GCs residing in it.

Observations reported here were obtained at the MMT Observatory, a joint facility of the University of Arizona and the Smithsonian Institution. MMT telescope time was granted, in part, by NOAO, through the Telescope System Instrumentation Program (TSIP). TSIP is funded by NSF. E.W.P. thanks the staff of the MMT Observatory for their professional support over multiple observing runs. E.W.P. also thanks the Harvard-Smithsonian Center for Astrophysics and its Institute for Theory and Computation for hosting him during an unexpectedly long pandemic visit. We thank Daniel Fabricant for sharing his flux calibration routine for MMT/Hectospec data. We also thank the anonymous referee for useful comments.

This paper is based on observations obtained with Mega-Prime/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 Centre 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.

C.L. acknowledges support from the National Natural Science Foundation of China (NSFC, grant Nos. 12173025, 11673017, 11833005, 11933003), 111 project (No. B20019), and Key Laboratory for Particle Physics, Astrophysics and Cosmology, Ministry of Education. A.L. is supported by Fondazione Cariplo, grant No. 2018-2329. A.L., T.H.P., and P.-A.D. acknowledge support from grant ANR-19-CE31-0022 (POPSYCLE) of the French Agence Nationale de la Recherche. L.V.S acknowledges support from the NASA ATP 80NSSC20K0566 and NSF CAREER 1945310 grants. A.J. and S.E. acknowledge support from ANID—Millennium Science Initiative—ICN12_009. M.G.L. was supported by the National Research Foundation grant funded by the Korean Government (NRF-2019R1A2C2084019). H.S.H. acknowledges the support from the Sejong Science Fellowship Program through the National Research Foundation of Korea (NRF-2021R1C1C2006790).

ORCID iDs

Youkyung Ko https://orcid.org/0000-0001-6333-599X

Eric W. Peng https://orcid.org/0000-0002-2073-2781

Patrick Côté https://orcid.org/0000-0003-1184-8114

Laura Ferrarese https://orcid.org/0000-0002-8224-1128

Chengze Liu https://orcid.org/0000-0002-4718-3428

Ariane Lançon https://orcid.org/0000-0002-7214-8296

Roberto P. Muñoz https://orcid.org/0000-0003-1743-0456

Thomas H. Puzia https://orcid.org/0000-0003-0350-7061

Karla A. Alamo-Martínez https://orcid.org/0000-0002-5897-7813

Laura V. Sales https://orcid.org/0000-0002-3790-720X
Felipe Ramos-Almendares https://orcid.org/0000-0002-2220-3736

Mario G. Abadi https://orcid.org/0000-0003-3055-6678

Myung Gyoon Lee https://orcid.org/0000-0003-2713-6744

Ho Seong Hwang https://orcid.org/0000-0003-3428-7612

Nelson Caldwell https://orcid.org/0000-0003-2352-3202

John P. Blakeslee https://orcid.org/0000-0002-5213-3548

Alessandro Boselli https://orcid.org/0000-0002-9795-6433

Jean-Charles Guillaume https://orcid.org/0000-0002-3263-8645

Pierre-Alain Duc https://orcid.org/0000-0003-3343-6284

Puragra Guhathakurta https://orcid.org/0000-0001-8867-4234

Stephen Gwyn https://orcid.org/0000-0001-8221-8406

Andrés Jordán https://orcid.org/0000-0002-5389-3944

Sungsoon Lim https://orcid.org/0000-0002-5049-4390

Rubén Sánchez-Janssen https://orcid.org/0000-0003-4945-0056

Elisa Toloba https://orcid.org/0000-0001-6443-5570

## References

Ahn, C. P., Alexandroff, R., Allende Prieto, C., et al. 2014, ApJS, 211, 17

Amorisco, N. C. 2019, MNRAS, 482, 2978

Andrews, B. H., & Martini, P. 2013, ApJ, 765, 140

Annibali, F., Bressan, A., Rampazzo, R., et al. 2007, A&A, 463, 455

Arnaboldi, M., Freeman, K. C., Gerhard, O., et al. 1998, ApJ, 507, 759

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

Beasley, M. A., Sharples, R. M., Bridges, T. J., et al. 2000, MNRAS, 318, 1249

Beasley, M. A., Trujillo, I., Leaman, R., et al. 2018, Natur, 555, 483

Bernstein, G. M., Nichol, R. P., Tyson, J. A., et al. 1995, AJ, 110, 1507

Bílek, M., Duc, P.-A., Cuillandre, J.-C., et al. 2020, MNRAS, 498, 2138

Binggeli, B., Popescu, C. C., & Tammann, G. A. 1993, A&AS, 98, 275

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

Binggeli, B., Tammann, G. A., & Sandage, A. 1987, AJ, 94, 251

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

Boselli, A., Fossati, M., Ferrarese, L., et al. 2018, A&A, 614, A56

Bovy, J., Hogg, D. W., & Roweis, S. T. 2011, AnApS, 5, 1657

Burstein, D., Faber, S. M., Gaskell, C. M., & Krumm, N. 1984, ApJ, 287, 586

Buzzoni, A., Arnaboldi, M., & Corradi, R. L. M. 2006, MNRAS, 368, 877

Cabrera-Ziri, I., & Conroy, C. 2022, MNRAS, 511, 341

Caldwell, N., Schiavon, R., Morrison, H., et al. 2011, AJ, 141, 61

Cappellari, M., Emsellem, E., Krajnović, D., et al. 2011, MNRAS, 413, 813

Cardiel, N., Gorgas, J., Cenarro, J., & Gonzalez, J. J. 1998, A&AS, 127, 597

Caso, J. P., De Bórtoli, B. J., Ennis, A. I., et al. 2019, MNRAS, 488, 4504

Chen, Y., & Gnedin, O. Y. 2022, arXiv:2203.00599

Chilingarian, I. V., Cayatte, V., Durret, F., et al. 2009, A&A, 486, 85

Ciardullo, R. 2010, PASA, 27, 149

Ciardullo, R., Durrell, P. R., Laychak, M. B., et al. 2004, ApJ, 614, 167

Cohen, J. G., Blakeslee, J. P., & Côté, P. 2003, ApJ, 592, 866

Cohen, J. G., Blakeslee, J. P., & Ryzhov, A. 1998, ApJ, 496, 808

Contini, E., De Lucia, G., Villalobos, Á., et al. 2014, MNRAS, 437, 3787

Cooper, A. P., D'Souza, R., Kauffmann, G., et al. 2013, MNRAS, 434, 3348

Cortesi, A., Arnaboldi, M., Coccato, L., et al. 2013, A&A, 549, A115

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

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

Duc, P.-A., Cuillandre, J.-C., Karabal, E., et al. 2015, MNRAS, 446, 120

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

Elson, R. A. W., & Santiago, B. X. 1996, MNRAS, 278, 617

Fabricant, D., Fata, R., Roll, J., et al. 2005, PASP, 117, 1411

Fabricant, D. G., Kurtz, M. J., Geller, M. J., et al. 2008, PASP, 120, 1222

Feldmeier, J. J., Mihos, J. C., Morrison, H. L., et al. 2004, ApJ, 609, 617

Fensch, J., Duc, P.-A., Lim, S., et al. 2020, A&A, 644, A164

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

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

Ferrarese, L., Côté, P., MacArthur, L. A., et al. 2020, ApJ, 890, 128

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

Forbes, D. A., & Remus, R.-S. 2018, MNRAS, 479, 4760

Gallazzi, A., Charlot, S., Brinchmann, J., et al. 2005, MNRAS, 362, 41

Gavazzi, G., Boselli, A., Donati, A., et al. 2003, A&A, 400, 451

Gonzalez, A. H., Zabludoff, A. I., Zaritsky, D., et al. 2000, ApJ, 536, 561

Graves, G. J., & Schiavon, R. P. 2008, ApJS, 177, 446

Greene, J. E., Janish, R., Ma, C.-P., et al. 2015, ApJ, 807, 11

Gu, M., Conroy, C., & Brammer, G. 2018, ApJL, 862, L18

Gu, M., Conroy, C., Law, D., et al. 2018, ApJ, 859, 37

Harris, W. E. 2009, ApJ, 703, 939

Harris, W. E., Harris, G. L., & Hudson, M. J. 2015, ApJ, 806, 36

Hernquist, L. 1990, ApJ, 356, 359

Hudson, M. J., Gillis, B. R., Coupon, J., et al. 2015, MNRAS, 447, 298

Hui, X., Ford, H. C., Ciardullo, R., et al. 1993, ApJ, 414, 463

Hwang, H. S., Lee, M. G., Park, H. S., et al. 2008, ApJ, 674, 869

Iodice, E., Sarzi, M., Bittner, A., et al. 2019, A&A, 627, A136

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

Johansson, P. H., Naab, T., & Ostriker, J. P. 2012, ApJ, 754, 115

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

Jordán, A., Peng, E. W., Blakeslee, J. P., et al. 2009, ApJS, 180, 54

Kartha, S. S., Forbes, D. A., Alabi, A. B., et al. 2016, MNRAS, 458, 105

Kenney, J. D. P., Tal, T., Crowl, H. H., et al. 2008, ApJL, 687, L69

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

Kirby, E. N., Cohen, J. G., Guhathakurta, P., et al. 2013, ApJ, 779, 102

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

Ko, Y., Lee, M. G., Park, H. S., et al. 2018, ApJ, 859, 108

Ko, Y., Lee, M. G., Park, H. S., et al. 2019, ApJ, 872, 202

Ko, Y., Lee, M. G., Park, H. S., et al. 2020, ApJ, 903, 110

Koleva, M., Prugniel, P., Ocvirk, P., et al. 2008, MNRAS, 385, 1998

Krick, J. E., & Bernstein, R. A. 2007, AJ, 134, 466

Kruijssen, J. M. D., Pfeffer, J. L., Crain, R. A., et al. 2019, MNRAS, 486, 3134

Kundu, A., Whitmore, B. C., Sparks, W. B., et al. 1999, ApJ, 513, 733

Kurtz, M. J., & Mink, D. J. 1998, PASP, 110, 934

Lambert, R. A., Rhode, K. L., & Vesperini, E. 2020, ApJ, 900, 45

Lamers, H. J. G. L. M., Kruijssen, J. M. D., Bastian, N., et al. 2017, A&A, 606, A85

Larsen, S. S., Brodie, J. P., Huchra, J. P., et al. 2001, AJ, 121, 2974

Lee, H.-c., Yoon, S.-J., & Lee, Y.-W. 2000, AJ, 120, 998

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

Lee, M. G., Park, H. S., Hwang, H. S., et al. 2010b, ApJ, 709, 1083

Lequeux, J., Peimbert, M., Rayo, J. F., et al. 1979, A&A, 80, 155

Li, C., Zhu, L., Long, R. J., et al. 2020, MNRAS, 492, 2775

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

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

Liu, C., Peng, E. W., Jordán, A., et al. 2011, ApJ, 728, 116

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

Longobardi, A., Arnaboldi, M., Gerhard, O., et al. 2015a, A&A, 579, A135

Longobardi, A., Arnaboldi, M., Gerhard, O., et al. 2015b, A&A, 579, L3

Longobardi, A., Arnaboldi, M., Gerhard, O., et al. 2018a, A&A, 620, A111

Longobardi, A., Peng, E. W., Côté, P., et al. 2018b, ApJ, 864, 36

Mannucci, F., Cresci, G., Maiolino, R., et al. 2010, MNRAS, 408, 2115

Martín-Navarro, I., La Barbera, F., Vazdekis, A., et al. 2015, MNRAS, 451, 1081

McDermid, R. M., Alatalo, K., Blitz, L., et al. 2015, MNRAS, 448, 3484

McElroy, D. B. 1995, ApJS, 100, 105

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

Mihos, J. C., Harding, P., Feldmeier, J., et al. 2005, ApJL, 631, L41

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

Mihos, J. C., Harding, P., Rudick, C. S., et al. 2013, ApJL, 764, L20

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

Mouhcine, M., Ferguson, H. C., Rich, R. M., et al. 2005a, ApJ, 633, 821

Mouhcine, M., Rejkuba, M., & Ibata, R. 2007, MNRAS, 381, 873

Mouhcine, M., Rich, R. M., Ferguson, H. C., et al. 2005b, ApJ, 633, 828

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

Murante, G., Giovalli, M., Gerhard, O., et al. 2007, MNRAS, 377, 2

Muratov, A. L., & Gnedin, O. Y. 2010, ApJ, 718, 1266

NASA/IPAC Extragalactic Database (NED) 2019, NASA/IPAC Extragalactic Database (NED), IPAC doi:10.26132/NED1

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

Oosterloo, T., & van Gorkom, J. 2005, A&A, 437, L19

Oser, L., Ostriker, J. P., Naab, T., et al. 2010, ApJ, 725, 2312

Park, H. S., & Lee, M. G. 2013, ApJL, 773, L27

Pastorello, N., Forbes, D. A., Usher, C., et al. 2015, MNRAS, 451, 2625

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

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

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, ApJS, 681, 197

Pfeffer, J., Kruijssen, J. M. D., Crain, R. A., et al. 2018, MNRAS, 475, 4309

Pota, V., Forbes, D. A., Romanowsky, A. J., et al. 2013, MNRAS, 428, 389
Powalka, M., Puzia, T. H., Lançon, A., et al. 2018, ApJ, 856, 84

Proctor, R. N., Forbes, D. A., & Beasley, M. A. 2004, MNRAS, 355, 1327

Purcell, C. W., Bullock, J. S., & Zentner, A. R. 2007, ApJ, 666, 20

Park, H. S., Lee, M. G., Hwang, H. S., et al. 2012, ApJ, 759, 116

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

Ramos-Almendares, F., Sales, L. V., Abadi, M. G., et al. 2020, MNRAS, 493, 5357

Recio-Blanco, A., de Laverny, P., Kordopatis, G., et al. 2014, A&A, 567, A5

Reina-Campos, M., Kruijssen, J. M. D., Pfeffer, J. L., et al. 2019, MNRAS, 486, 5838

Reina-Campos, M., Trujillo-Gomez, S., Deason, A. J., et al. 2021, arXiv:2106.07652

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

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

Sánchez-Janssen, R., Côté, P., Ferrarese, L., et al. 2019, ApJ, 878, 18

Schiavon, R. P., Rose, J. A., Courteau, S., et al. 2004, ApJL, 608, L33

Schiavon, R. P., Rose, J. A., Courteau, S., et al. 2005, ApJS, 160, 163

Scodeggio, M., Giovanelli, R., & Haynes, M. P. 1998, AJ, 116, 2738

Searle, L., & Zinn, R. 1978, ApJ, 225, 35

Smith, R. J., Lucey, J. R., Hudson, M. J., et al. 2009, MNRAS, 392, 1265

Spolaor, M., Kobayashi, C., Forbes, D. A., et al. 2010, MNRAS, 408, 272

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

Tamura, N., Sharples, R. M., Arimoto, N., et al. 2006, MNRAS, 373, 601

Thomas, D., Maraston, C., & Bender, R. 2003, MNRAS, 339, 897

Thomas, D., Maraston, C., Bender, R., et al. 2005, ApJ, 621, 673

Thomas, D., Maraston, C., & Johansson, J. 2011, MNRAS, 412, 2183

Trager, S. C., Worthey, G., Faber, S. M., et al. 1998, ApJS, 116, 1

Tremonti, C. A., Heckman, T. M., Kauffmann, G., et al. 2004, ApJ, 613, 898

Trujillo-Gomez, S., Kruijssen, J. M. D., Reina-Campos, M., et al. 2021, MNRAS, 503, 31

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

Villaume, A., Foreman-Mackey, D., Romanowsky, A. J., et al. 2020, ApJ, 900, 95

Vogelsberger, M., Genel, S., Springel, V., et al. 2014a, MNRAS, 444, 1518

Vogelsberger, M., Genel, S., Springel, V., et al. 2014b, Natur, 509, 177

Wang, Q., Peng, E. W., Blakeslee, J. P., et al. 2013, ApJ, 769, 145

Whitmore, B. C., Sparks, W. B., Lucas, R. A., et al. 1995, ApJL, 454, L73

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

Williams, B. F., Dalcanton, J. J., Bell, E. F., et al. 2012, ApJ, 759, 46

Worthey, G., Faber, S. M., Gonzalez, J. J., & Burstein, D. 1994, ApJS, 94, 687

Worthey, G., & Ottaviani, D. L. 1997, ApJS, 111, 377

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

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