Nonlinear Color-Metallicity Relations of Globular Clusters. X. Subaru/FOCAS Multi-object Spectroscopy of M87 Globular Clusters

We obtained spectra of some 140 globular clusters (GCs) associated with the Virgo central cD galaxy M87 with the Subaru/FOCAS MOS mode. The fundamental properties of GCs such as age, metallicity and $\alpha$-element abundance are investigated by using simple stellar population models. It is confirmed that the majority of M87 GCs are as old as, more metal-rich than, and more enhanced in $\alpha$-elements than the Milky Way GCs. Our high-quality, homogeneous dataset enables us to test the theoretical prediction of inflected color$-$metallicity relations (CMRs). The nonlinear-CMR hypothesis entails an alternative explanation for the widely observed GC color bimodality, in which even a unimodal metallicity spread yields a bimodal color distribution by virtue of nonlinear metallicity-to-color conversion. The newly derived CMRs of old, high-signal-to-noise-ratio GCs in M87 (the $V-I$ CMR of 83 GCs and the $M-T2$ CMR of 78 GCs) corroborate the presence of the significant inflection. Furthermore, from a combined catalog with the previous study on M87 GC spectroscopy, we find that a total of 185 old GCs exhibit a broad, unimodal metallicity distribution. The results corroborate the nonlinear-CMR interpretation of the GC color bimodality, shedding further light on theories of galaxy formation.


INTRODUCTION
Globular clusters (GCs) are an old simple stellar population (SSP) with a small internal dispersion in age and chemical composition, and provide powerful clues to help understand the formation histories of their host galaxies. The exciting discovery of GC color bimodality in many luminous early-type galaxies (e.g., Geisler et al. 1996;Kundu & Whitmore 2001;Larsen et al. 2001;Peng et al. 2004;Harris et al. 2006;Peng et al. 2006;Lee et al. 2008;Jordán et al. 2009;Sinnott et al. 2010;Liu et al. 2011;Blakeslee et al. 2012;Brodie et al. 2012;Pota et al. 2013;Kartha et al. 2014;Richtler et al. 2015;Cho et al. 2016;Kartha et al. 2016;Caso et al. 2017;Harris et al. 2017;Choksi et al. intrinsic metallicity and its proxy, colors, holds the key to the color bimodality phenomenon. Paper I demonstrated that conversions to colors from the unimodal metallicity distribution of old GCs (> 10 Gyr) can be naturally achieved through nonlinear CMRs. Nonlinear relations between [Fe/H] and colors are also reported observationally (Peng et al. 2006;Blakeslee et al. 2012;Chies-Santos et al. 2012;Cho et al. 2016).
The nonlinear-CMR scenario was scrutinized in a subsequent series of papers (Yoon et al. 2011a(Yoon et al. ,b, 2013Kim et al. 2013;Chung et al. 2016;Kim & Yoon 2017;Lee et al. 2019;Lee et al. 2020, hereafter Papers II, III, IV, V, VI, VII, VIII, and IX). In Papers II and IV, the shapes of CMRs according to different colors were examined for the GC systems of M87 and M84, respectively. They showed that the u-band-related colors in particular, are useful diagnostics of the underlying metallicity distribution function. Different combinations of ugz bands from Hubble Space Telescope (HST) multiband photometry yielded histograms of varying morphologies, which is consistent with the nonlinear-CMR assumption. Paper III focused on the analysis of GC metallicities and their link to the host galaxy. It was shown that when transformed through nonlinear CMRs, GCs and halo field stars share similar metallicity distribution shapes. This suggests a connected origin of GCs and halo field stars, contrary to previous views. Papers V and VII examined the spectra of old, spheroidal GCs in the M31 and NGC 5128, respectively. The observed distributions of high-quality spectroscopic absorption-line indices of M31 GCs (Caldwell et al. 2011) and NGC 5128 GCs (Beasley et al. 2008;Woodley et al. 2010) were reproduced successfully using the nonlinear metallicityto-index transformation, exactly analogous to the view that the nonlinear metallicity-to-color transformation is responsible for the photometric color bimodality of GCs. A further theoretical spectroscopic study, Paper VI, showed that the predictions from the Yonsei Evolutionary Population Synthesis (YEPS) models (Chung et al. 2013;Chung et al. 2017) exhibit nonlinearity in the Ca II triplet (CaT ; 8498, 8542, and 8662 A) index versus metallicity plane. The CaT is a wellknown metallicity indicator for stellar populations and is used to derive metallicities of extragalactic GCs assuming a linear relation between CaT and metallicity (e.g., Foster et al. 2010;Brodie et al. 2012;Usher et al. 2015). They also showed that the conversion via nonlinear CaT −metallicity from a unimodal metallicity distribution reproduces the observed CaT distribution well. More recently, two theoretical photometric studies also reinforced the nonlinearity scenario for the GC color bimodality. Paper VIII showed that the individual color distributions of 78 GC systems in early-type galaxies from the ACS Virgo and Fornax Cluster Surveys are correctly reproduced based on the nonlinear-CMR scenario. Paper IX showed that the difference between blue and red GCs in M49, M60, M87, and NGC 1399 in terms of the radial surface number density profile arises naturally from the nonlinear metallicity-to-color conversion together with the observed radial metallicity gradient of the GC systems.
This paper is the 10th in the nonlinear-CMR series and adopts the same line of analysis in exploring and determining more precise forms of GC CMR. We use a homogeneous set of high-quality M87 GC spectra with wide metallicity ranges, obtained with the Subaru/FOCAS MOS mode. M87 is a cD galaxy at the center of the Virgo cluster harboring a population of ∼ 10 4 GCs. The previous spectroscopic study of its GC system (Cohen et al. 1998, hereafter C98) reported ages and metallicities for ∼ 100 GCs. Our sample GCs are located similarly within 10 ′ of galactic radius and have nine GCs in common with the C98 catalog. C98 spectroscopic GCs are brighter than ours by 0.5−1.0 mag in the V band. Recently, Villaume et al. (2019) reported a new metallicity measurement on M87 GCs via full-spectrum SSP model fitting using the C98 GC in which they obtained systematically higher metallicities compared to the previous findings on M87 (see the discussion).
The paper is outlined as follows. Section 2 presents the observation and reduction of the data used in the study. Section 3 describes the measurements of radial velocities and selection of member GCs. Section 4 provides the measurement of Lick/IDS absorption-line indices (Worthey 1994;Worthey & Ottaviani 1997, hereafter W94 and W97, respectively), followed by the inferred age, metallicity, and α-element abundance of the M87 GCs in Section 5. In Section 6, we present the derived CMRs of the clusters. The implications of our findings are discussed in Section 7.

OBSERVATIONS AND DATA REDUCTION
The spectroscopic GC targets are selected from the photometric data on the M87 GC candidates obtained with Suprime-Cam on the 8.2 m Subaru telescope (Tamura et al. 2006a;Tamura et al. 2006b). In Figure  1, the V versus V − I color−magnitude diagram (CMD; upper panel) and the color histogram (lower panel) show that the selected clusters cover the typical color range of GCs (0.8 < V − I <1.4) and exhibit a clear color bimodality. Figure 2 shows the locations of the GCs as blue (V − I < 1.1) and red (V − I > 1.1) filled circles. In the inner target fields (< 10 ′ ) where GC density is very high, we prioritized the red clus-ters with V = 21.5 − 22.0 mag because, at larger distances, color bimodality is less clear due to the decreasing contribution of red GCs. The magnitude range is chosen in order to avoid the possible effect of the known blue tilt of the M87 GC system at V ≤ 21.5 (Strader et al. 2006). The phenomenon of a lack of bright blue clusters is commonly observed in massive elliptical galaxies (Ostrov et al. 1993;Dirsch et al. 2003;Harris et al. 2006;Mieske et al. 2006;Strader et al. 2006;Harris et al. 2009;Mieske et al. 2010).
We performed spectroscopy on a total of 157 M87 GC candidates with the Multi-Object Spectrography (MOS) mode of the Faint Object Camera and Spectrograph (FOCAS) on the Subaru telescope on 2007 March 16 − 19. We used a configuration of a 300B grism and L600 filter resulting in a dispersion of 1.35Å/pixel and a spectral coverage of 3700 − 6000Å. Object spectra were obtained through 0. ′′ 5 slits, rendering an overall spectral resolution of R ∼ 800 or 5.4Å, which is calculated using the spectral dispersion (1.35Å/pixel) and the FWHM of a sky line 4.0 pixel. Seeing was on average 1 ′′ , ranging between 0. ′′ 9 and 1. ′′ 5 over the four nights. FOCAS has a circular 6 ′ field of view and the four masks contained between 28 and 46 objects ( Table 1). The typical exposure time for the M87 fields was 1800 s, which adds up to the total integration time of ∼ 22 hr (Table 1). In addition, we observed a set of calibration stars including two flux standard stars (Feige34 and HZ44), three RV standard stars (HD68874, HD69148, and HD160952), and 26 Lick standard stars with a long-slit mode. Dome flats and arc frames for both long-slit and MOS observations were taken in between exposures for object spectra.
We used basic spectroscopic data reduction procedures (bias subtraction, cosmic-ray removal, flatfielding, field distortion correction, flexure effect correction, wavelength calibration, sky subtraction, spectrum extraction, extinction correction, and flux calibration) using the FOCASRED package. We applied 3 pixel on-chip binning in the spatial direction and normalized the spectra by the median flux density at 4900 − 5100 A. Finally, to further explore the CMRs of M87 GCs, we have also cross-identified our spectroscopic targets with the M − T 2 Washington photometry data. We obtained M and T 2 images of M87 GCs with the Mosaic II CCD imager of the 4 m Blanco telescope at the Cerro Tololo Inter-American Observatory (CTIO). The Mosaic II CCD consists of eight 2048 × 4096 array with a pixel scale of 0.27, yielding a field of view of 36 ′ × 36 ′ . Images were preprocessed through IRAF 1 using the MSCRED package (Valdes 1998). The standard data reduction was then performed, and the photometry of the objects was obtained using the DAOPHOT II/ALLSTAR package (Stetson 1987). We used standard stars by Landolt (1992) and Stetson (2000) for the photometric calibration. Eight standard fields, each containing 15-25 standard stars with a wide color range, were taken at various airmasses. Reddening correction and astrometry were done using the reddening maps by (Schlegel et al. 1998) and the USNO-B1.0 catalog (Monet et al. 2003), respectively.

RADIAL VELOCITY MEASUREMENT
Among our selected sample cluster candidates, we now identify bona fide cluster members through the measurement of RVs of the candidate GCs. RVs of the 157 GC candidates are measured by the Fourier cross-correlation technique (Tonry & Davis 1979), using the FXCOR task in the IRAF RV package. We fit the continuum of the spectra using the spline fit with 3σ clipping, then subtract the fit from the original spectra of GC candidates. For the reference spectra required for this method, we observed three velocity standard stars of spectral types G8III during the observing run. Then the continuumsubtracted spectra of GC candidates are cross-correlated with those of RV stars. Each velocity standard star is separately used, yielding three independent values of RVs for GCs. The derived velocity values using three different standard stars are consistent with each other within 1σ, and much smaller than the individual measurement errors. Because the cross-correlation technique is rather sensitive to the noisy part of an input spectrum, which can smear the correlation signal, the choice of wavelength range to be cross-correlated is important. We choose the wavelength region of 4000−5500Å for the task in order to avoid the relatively low signal-to-noiseratio (S/N) region of our object spectra below ∼ 4000 A, and a strong sky emission-line region around 5600Å. This range includes 15 − 18 Lick features. In obtaining the final velocity values for GCs, the cross-correlated range of wavelengths was adjusted on an individual object basis. Given the consistency between the velocity values derived with different velocity standard stars, the final RVs of GCs are measured using the RV standard star HD 68874 based on the quality of the standard star data. The uncertainties of the RVs vary with the mask fields due to the S/N variation of the fields. The uncertainties range mostly from 35 − 73 km s −1 , with the typical error of ∼ 61 km s −1 within 10 ′ fields.
The typical velocities of M87 GCs are known to be 200 − 2550 km s −1 based on the previous kinematic studies of the system (C98; Hanes et al. 2001). Because objects with RV < 200 km s −1 are almost definitely foreground Galactic stars, we consider the velocity criterion of 200 < RV < 2550 km s −1 for the GCs associated with M87. Except for obvious contamination (emission-line galaxies and/or AGN) and objects whose spectra are useless due to bad columns, we find that 130 GCs satisfy this criterion. The measured velocities of all cluster candidates are presented in Table 2 along with the corresponding magnitude and color information by Tamura et al. (2006a); Tamura et al. (2006b) in the fourth and fifth columns.
The histogram of RVs of the confirmed GCs is shown in Figure 3. The RV distribution of M87 GCs has double peaks, and the clusters are distributed roughly symmetrically around the mean velocity of M87 GCs, 1299 km s −1 . This value is in good agreement with the known RV of M87 (1307 km s −1 , C98; Hanes et al. 2001;Strader et al. 2011, hereafter S11) at a distance of 16.1 Mpc.
Some of the GCs have previous measurements by C98 and S11. Figure 4 illustrates these GCs in comparisons between each set of measurements. The left panel shows nine GCs in common with the C98 catalog, and the right panel shows 39 GCs in common with the new velocity measurements by S11. The velocities of these common GCs are in general consistent. The dotted line in each panel represents the least-squares fit to the data, and the rms of both fits are found to be comparable to the measurement errors. One exception is F233 (RV = 962 km s −1 ) which has measurements in all three catalogs. The velocity value for this cluster by C98 is 1762 km s −1 . The cause of this discrepancy is unknown. Our value for this cluster is in good agreement with the more recent measurement by S11.
The Lick/IDS system of absorption-line indices provides a reference frame in which one can derive metallicities, ages, and abundance ratios for stellar populations. The Lick/IDS system was originally created by Burstein et al. (1984), and has been continuously revised and improved over the years (W94; Trager et al. 1998). The system has also been tried and developed in SSP modeling, with their variations with stellar age, metallicity, and α-element abundance considered (e.g., Bruzual & Charlot 2003;Thomas et al. 2003;Chung et al. 2013). The observed Lick index values are compared to the SSP Lick indices to estimate ages, metallicities, and [α/Fe] ratios. Figure 5 shows three representative spectra of GCs with RV = 0 km s −1 and the marked locations of the absorption features (gray solid lines). The metallicity [Fe/H] of the spectra increases from the top to bottom panel. Following the wavelength-dependent Lick spectral resolution, we smoothed our spectra (red solid lines) to match the original Lick resolutions (W97). All indices and errors were measured using the C ++ program Indexf (Cardiel et al. 1998). In this study, we adopt the passband and pseudo-continuum definitions by W94 and W97. The object spectra and error spectra were provided as inputs for Indexf Lick index measurements. Indexf generates Lick index uncertainties by performing 100 Monte Carlo simulations; indices are measured on the Poisson noise added spectra, and the 1σ standard deviation of the measured index is taken as the Lick index uncertainty.
In Figure 6, we calibrate our spectra to the Lick system by comparing our measured indices for the 26 Lick standard stars to the Lick/IDS measurements. The offsets between the two measurements were then applied to the Lick indices of our 130 M87 GCs listed in Table  2. The measured Lick indices and their uncertainties are given in Tables 3 and 4, respectively. In addition to the above Lick indices, we have also combined individual indices to create a more robust composite index, such as (Thomas et al. 2003) and We show the comparisons of our Lick index measurements with those by C98 in Figure 7. Among the 150 GCs obtained by C98, 10 GCs are also on our slit masks and the Lick indices were successfully measured for 9 of these. After excluding one GC with inconsistent RV values, Lick indices of eight common GCs are compared. While Mgb shows good agreement between the two data sets, the discrepancies are apparent for Fe and Balmer indices. We note that the Lick indices by C98 are measured using the older definitions by Burstein et al. (1984), Faber et al. (1985) and Gorgas et al. (1993), while ours adopts the definitions by Trager et al. (1998). This may be partly responsible for some of the larger rms values.

AGE, METALLICITY, AND ELEMENTAL ABUNDANCE
In this section, we explore the stellar population properties-ages, metallicities, and elemental abundance ratios-for our sample GCs via various methods and check for consistency between the results. In Section 5.1, we first examine the general ranges and trends of GC age, metallicity, and abundances in the index−index diagnostic diagrams by comparing the observations to the SSP model grids. Then, in Section 5.2, we obtain the purely empirical metallicity values, which are independent of models, using the most recent comprehensive catalog of Milky Way globular clusters (MWGCs, Kim et al. 2016). These MWGCs are old and of solar abundance ratios. Next, in Section 5.3, in order to assign ages, metallicities, and elemental abundance ratio values to each individual GC we adopt the method of multi-index fitting to a set of SSP models via χ 2 minimization described by Proctor et al. (2004). The derived values in Section 5.3 are used in the analysis for the rest of the paper.

Absorption-line Diagnostic Plots
The index−index diagnostics provide a simple, qualitative picture on the distributions of ages, metallicities, and elemental abundance ratios of a GC system. In Figure 8, we present plots with selected indices as a function of the metallicity-sensitive composite index [MgFe] ′ (Equations 1). The combination of indices containing the three Balmer indices Hβ, Hγ A , and Hγ F represent the age−metallicity distributions of the clusters. The superimposed grids are the model predictions from the Yonsei Evolutionary Population Synthesis (YEPS) model (Chung et al. 2013)  To estimate the abundance ratios [α/Fe] of the clusters, we choose tracers of α-elements (Mgb) and ironpeaked elements ( Fe , Equations 2) (lower-right panel). α-elements and more than two-thirds of iron-peak ele-ments are produced by Type II and Type Ia supernovae (SNe), respectively. By virtue of the distinct characteristics of each SN type in terms of timing and duration of events, the ratio of [α/Fe] provides clues to the star-formation history of the host. Enhanced [α/Fe] (∼ 0.3) is commonly found for GCs in early-type galaxies (Worthey et al. 1992;Matteucci 1994;Trager et al. 2000;Thomas et al. 2005), indicating rapid star formation. The observational points are denoted with the symbols as before, and the 12 Gyr model GCs with varying [α/Fe] ratios of 0.0 (dotted line), 0.3 (solid line), and 0.6 (dashed line) are overlaid. The majority of M87 GCs lie within the range of [α/Fe] = 0.0−0.3 which is consistent with the previous findings on the extragalactic GC system for early-type galaxies. At the low-metallicity end, it is difficult to draw a conclusion on the distribution of [α/Fe] due to convergence of the model grids as well as somewhat large observational uncertainties.
We next examine the abundance pattern of carbon and nitrogen through a set of index−index diagnostics of M87 GCs using CN 2 and G4300 indices ( Figure 9). In the left panel, the CN 2 distribution of M87 GCs are shown as light-gray filled squares, and the clusters with higher S/N are plotted as dark gray squares. Although there is a large scatter of M87 data compared to the MWGCs (open squares; Kim et al. 2016), a general trend with metallicity is obvious. The overplotted model grids represent the varying [α/Fe] ratios of solar (gray) and 0.3 dex (black), respectively. CN 2 exhibits systematically stronger index strengths for both MWGCs and M87 GCs across the full metallicity range with respect to the model predictions. This indicates carbon and/or nitrogen overabundances for the two GC systems, the trend of which seems to be common in typical GC systems (Puzia et al. 2005;Cenarro et al. 2007;Beasley et al. 2008;Woodley et al. 2010). In the right panel, G4300 is known to be mainly sensitive to carbon (Tripicco & Bell 1995) and shows a smaller offset between model predictions and the observation.
The CN-strong features observed in GCs are ascribed to several possible sources, including pollution of gas by intermediate-mass AGB stars (Cottrell & Da Costa 1981;D'Antona & Ventura 2007), early enrichment achieved by the stellar wind of fast-rotating massive stars (Maeder & Meynet 2006;Prantzos & Charbonnel 2006;Decressinet al. 2007), and difference in surface gravity (Mucciarelli et al. 2015;Lee 2016;Gerber et al. 2018), opacity (MacLean et al. 2018) and stellar mass function (Renzini 2008;Kim & Lee 2018). The nitrogen overabundance may also be attributed to the two generations of the constituent stars of GCs (Carretta et al. 2008;D'Ercole et al. 2008). In this scenario, second-generation stars are formed from the material preenriched by supernovae explosions. Because the observed G4300 indices with higher S/N do not seem to stray far from the models, we can assume that carbon enhancement is not likely. Hence, overenrichment of nitrogen in M87 GCs may be a source of the significant discrepancy between the observation and model. Due to large uncertainties, we cannot conclusively determine the abundance pattern of M87 GCs, but the overall trends of CN 2 and G4300 against the metal-sensitive index [MgFe] ′ indicate an overabundance in carbon and nitrogen.

Empirical Metallicities
We employ a simple, model-independent approach to assigning metallicities to individual M87 GCs in which MWGCs are used as metallicity calibrators. It is noted in Beasley et al. (2008) and Strader & Brodie (2004) that, although mainly tied to the Zinn & West (1984) metallicity scale, there is no metallicity scale universally agreed upon. In this study, we refer to our empirically derived metallicity as [M/H] which may not strictly reflect either [Fe/H] or overall metallicity ([Z/H] (see Beasley et al. (2008) for a more detailed discussion). We use spectroscopic data on 53 MWGCs (Kim et al. 2016) which are calibrated onto the Lick system and thus directly comparable to our M87 GC data.
We adopt the method that closely follows that of Beasley et al. (2008). In order to derive empirical metallicities, we obtain correlations between metallicities provided by the Harris catalog (Harris 1996) and the measured Lick indices of 53 MWGCs (Kim et al. 2016). To better define the relationships that exhibit a hint of nonlinearity, we fit the multiorder polynomial to the data using orthogonal distance regression (ODR;Jefferys et al. 1988). In most cases, the correlations are characterized by second-order polynomial, except for Hβ and CN indices, which are better defined by higher-order correlations. Figure 10 illustrates the resulting fits for the 18 Lick indices and the fitted coefficients are listed in Table  5. Nine indices are chosen for the derivation of empirical metallicities based on both the goodness of the ODR fits on the MWGCs and the index measurement qualities for M87 GCs. These indices include Fe4383, Fe4531, Hβ, Fe5015, Mgb, Mg 2 , Fe5270, Fe5335, and Fe5406. As shown previously in Figure 8, young GCs cannot be readily separated due to the possible effect of hot horizontal-branch stars present in the GCs. Taking this into consideration, to examine the properties of the old members of the M87 GC system exclusively, we leave out four clusters that are found to be ≤ 8 Gyr in all three Balmer−[MgFe] ′ planes ( Figure 8). As the spectroscopic ages of M87 GCs and stellar population are in general old ( > 12 Gyr, Cohen et al. (1998);Kuntschner et al. (2010)), we assume the rest are old GCs. For 83 old GCs with sufficient quality (10 < S/N < 40) for an abundance analysis (Figure 11), we obtain empirical metallicities [M/H] using the fitted coefficients for the chosen indices and list the values in Table 3.

Multi-index Fits via χ 2 Minimization
For a more quantitative analysis on the metallicities, ages, and [α/Fe] of M87 GCs, we perform χ 2 multiindex fitting with the YEPS SSP models (Chung et al. 2013). This technique, introduced and explained in detail by Proctor et al. (2004), simultaneously fits a set of Lick indices to the model grids of [Fe/H], ages, and [α/Fe] until χ 2 is minimized. The method has the advantage in that by using multiple absorption-line indices, it is more robust against individual index uncertainties. We adopt the same method of rejecting the deviant indices and recalculating the χ 2 statistics as described in Proctor et al. (2004). Note that we exclude several indices (e.g., CN 1 , CN 2 , Ca4227 and G4300) from the initial fitting process because of their inconsistency with respect to the models. We determined 13.5 Gyr to be the best-fit model age parameter to represent the old clusters in M87 and thus fixed the model as such in our multi-index fitting process to obtain the metallicities [Fe/H] and abundance ratios [α/Fe]. We perform a χ 2 routine with the same 83 GCs as the previous section. The [Fe/H] values are compared to the empirical metallicities derived in Section 5.2 ( Figure 12) and found to be in good agreement with each other. Because the empirical [M/H] is set by MWGCs, we might expect offsets for GCs with [α/Fe] higher than MWGCs' abundance ratios (see the bottom-right panel of Figure 8). GCs with higher α-element abundance ratios ([α/Fe] > 0.5) are denoted with gray squares and they exhibit larger inconsistencies in the metal-poor region. These GCs are excluded from the least-squares fits to the data. The derived values of ages, metallicities, and abundance ratios for 87 GCs in this section are used in the analysis for the rest of the paper.
6. TESTS FOR NONLINEAR-CMR SCENARIO FOR GC COLOR BIMODALITY

Color−Metallicity Relations
Our high-quality, homogeneous dataset enables us to test the theoretical prediction of inflected CMRs. This figure is similar to Figure 2 of Paper II where the model CMR was explored for a combined heterogeneous data set of GCs in M87, M49, and the MW, but it now uses the by far more homogeneous dataset of M87 GCs only. We find that the CMRs of M87 GCs in both metallicity−(V − I) and metallicity−(M − T 2) planes display wavy features where slopes of the relations change at around [Fe/H] = −1. The observed significant inflection along the optical CMRs supports the notion that the nonlinear CMRs between intrinsic metallicity and its proxy, colors, holds the key to the color bimodality phenomenon. The effect of this nonlinear conversion from metallicities to colors will be demonstrated in Section 6.3.

Metallicity Distribution Functions
Here we combine our metallicities derived with the SSP model-fitting method with C98 metallicities in order to have a better understanding of the intrinsic metallicity distribution function (MDF) shape of the M87 GC system. The left panel of Figure 14 illustrates the transformation of C98 metallicities to our M87 metallicities. The C98 catalog comprises bright GCs, which occupy regions similar to as our sample within 10 ′ of the galactic radius ( Figure 2). The solid line represents the orthogonal least-squares fit, and the derived transformation is given by In the right panel, the distribution of the M87 GC metallicities from this study is shown as a red hatched histogram. The blue hatched histogram displays the resulting metallicity distribution of C98 GCs with high-S/N C98 shifted to the same scale as that of ours. The gray filled histogram represents the combined metallicities of Subaru and Keck spectroscopy. The combined list of M87 GCs contains 98 higher-S/N GCs from C98 and 87 GCs from our study. For the GCs present in both catalogs, we take the values from our study. This is one of the largest spectroscopic data on M87 GCs and of quality that allows reliable abundance study on the system. A Gaussian mixture model (GMM) test on the observed MDFs (Table 7) gives a high probability of the distribution being a unimodal Gaussian one (P (χ 2 ) = 0.246).

Color Distribution Functions
In Figure 15 we simulate the M87 GC colors and compare them with the observations. The top row presents the observations and the best-fit model predictions (13.5 Gyr and [α/Fe] = 0.3) for the CMRs of M87 GCs. The models trace the features present in the observations well. In order to demonstrate the nonlinear-CMR hypothesis at work, we carry out the conversion process from metallicities to colors via our theoretical metallicity−color relations. For an MDF, we make a simple assumption of a single Gaussian distribution of 10 6 model GCs (top row, along the y-axes) to avoid small number statistics. The assumed mean [Fe/H] and σ([Fe/H]) of −0.85 and 0.5 are similar to those of observed distribution in Figure 13. In the second row, we plot the color−magnitude diagrams of 5000 randomly selected model GCs for the V − I and M − T 2 colors as an additional aid to visualize the simulated color distributions. We take observational uncertainties into account in the simulations. The inflection in our theoretical [Fe/H]−color relations has the effect of projecting the equidistant metallicity intervals onto wider color intervals, causing scarcity in the color domain near the quasi-inflection point on each [Fe/H]−color relation. As a result, the division into two groups of model colors are visible, which agrees with the observations ( Figure  1). The resultant color histograms (third row) show bimodal distributions with stronger blue peaks, and with the paucity at quasi-inflection points translated to the dip positions for both colors. Therefore, the theoretical prediction and the observed distributions of GC colors (bottom row) are consistent with each other. This result is also in good agreement with that of Paper II on the general behaviors of M87 GCs.
In order to quantify the visual impression of the distributions, we perform a GMM test on the simulated and observed color histograms. Table 7 presents the resulting GMM outputs. Bimodality is preferred for color distributions of the observed and simulated GCs (P (χ 2 ) = 0.001). The fraction of GCs assigned to the blue group is slightly higher in simulations. We add vertical dotted lines in Figure 15 at the blue and red peaks for the observed histograms (fourth row) to better guide the eye for a comparison with the simulation (third row). We attribute the slight offset in red peak locations for M − T 2 between the observed and simulated color histograms to several contributing factors, such as different choices of model ingredients involved. Hence, we put more weight on the consistent histogram morphology between the observations and simulations.

SUMMARY AND DISCUSSION
We have examined the fundamental stellar population properties of M87 GCs and reported spectroscopic metallicity estimates using the homogeneous spectro-scopic data. We have adopted empirical and theoretical methods in estimating the individual GC metallicity values and found them consistent with each other (Figure 12). The empirical metallicity was obtained using the most recent comprehensive catalog of MWGCs (Kim et al. 2016). The theoretical metallicity was based on multi-index fitting to a set of SSP models via χ 2 minimization described by Proctor et al. (2004).
The nonlinear-CMR scenario for GC color bimodality has been tested in the previous series of studies (Papers I ∼ IX) using multiple heterogeneous data sets of extragalactic GC systems, both photometric and spectroscopic, and found to be successful in explaining the bimodality phenomenon. Here we have used homogeneous spectroscopic data of the M87 GC system to continue to examine and test the nonlinear-CMR scenario. We have shown that M87 GCs exhibit nonlinearity in CMRs ( Figure 13). We have also shown that the MDF of the M87 GC system is characterized by a broad, unimodal distribution (Figure 14). In order to strengthen the statistical significance of the MDF, we have combined the metallicity measurements of our study and those of C98 with high S/N, making one of the largest catalogs of spectroscopically derived metallicity on M87 GCs. This combined distribution of metallicities derived with the spectroscopic data is consistent with the GC MDFs of several elliptical galaxies inferred from various colors through nonlinear CMRs (Papers II, III, and IV). Finally, the observed bimodal color distributions are reproduced using our nonlinear metallicity-to-color transformation scheme (Figure 15), bringing our results in good agreement with the findings of the previous studies.
Recently, Villaume et al. (2019) reported a new metallicity measurement on M87 GCs. Via full-spectrum SSP model fitting, they obtained a bimodal MDF with systematically higher metallicities compared to the previous findings on M87 (C98, Peng et al. (2006)). Their Figure 10 compares the M87 CMR with the literature (Peng et al. (2006)) and the MWGCs, and they found that the shapes of CMRs differ significantly. Villaume et al. suggested that the difference in the metal-poor region might be due to the age−metallicity degeneracy or the effect of blue horizontal-branch stars. Given that many extragalactic GCs are found to be old (e.g., C98; Forbes et al. 2001;Puzia et al. 2005) and that the YEPS models employ a realistic treatment of horizontalbranch stars in stellar population modeling, we do not consider them major causes for the discrepancy. We suspect the distinct CMRs might be the result of a different approach to deriving metallicities. Fahrion et al. (2020) compares GC metallicities of the galaxies in the Fornax cluster obtained with different techniques and various sets of constraints in their appendix. For example, the CMRs obtained by full-spectrum fitting and by using line-strength indices clearly exhibit offsets where GCs have different metallicities at fixed colors (see their Figure A.1).
In the conventional scenarios of GC formation, blue GCs come from the accretion of low-mass neighboring galaxies. However, whether the accretion of metal-poor GCs is solely responsible for blue GCs is still an open question. A piece of evidence against the accretion scenario is the observed sharp blue peaks of giant elliptical galaxies. The mean colors of GC systems correlate tightly with host galaxies' mass (Larsen et al. 2001;Strader & Brodie 2004;Peng et al. 2006). Such correlation suggests that accreted galaxies used to have a narrow mass range to form the sharp blue peak, which is unlikely. Further counterevidence is the significant fraction of blue GCs in galaxies in lower-density and isolated regions (Peng et al. 2006;Cho et al. 2012), where galaxies have few accompanying galaxies and thus it is difficult to acquire sufficient metal-poor GCs via accretion. Moreover, in our study, a unimodal shape of GC MDF is obtained in the inner main body of M87 (< 6.3 R eff ), where accreted GCs should be less populated than in the remote halo (see Paper III for further discussion).
Current hierarchical models of galaxy formation predict that the emergence of a massive galaxy involved a large number of protogalaxies, which may leave little room for the existence of merely two GC subpopulations in galaxies such as M87. In these models, the chemical evolution in the early stage of galaxy formation took place in a rather quasi-monolithic manner, which naturally produces a broad, unimodal MDF of GCs. We conclude from our results in this study, together with the evidence offered from our series of observational and theoretical studies (Papers I ∼ IX), that in the case of massive giant galaxies like M87, GCs are formed on a relatively short timescale via early, vigorous star formation (e.g., Paper III, Chiosi et al. 2014).     Beasley et al. 2008) and by using SSP models (x-axis; Proctor et al. 2004). The solid and dotted lines represent the orthogonal least-squares fits to the black squares and the one-to-one relations, respectively. Gray squares are the GCs with higher alpha-elemental abundance ratios ([α/Fe] > 0.5.).    Note-The full table is available in the online version.    Note-a The resulting significance of the GMM test expressed as a P-value. b The number ratios between the blue and red GC groups. c The mean and the standard deviations of the blue and red GC groups.  Note-a The resulting significance of the GMM test expressed as a P-value. b The number ratios between the blue and red GC groups. c The mean and the standard deviations of the blue and red GC groups.