The Assembly History of M87 through Radial Variations in Chemical Abundances of Its Field Star and Globular Cluster Populations

We model the spectra with an updated version of the absorption-line fitter (alf, Conroy & van Dokkum 2012; Choi et al. 2014; Conroy et al. 2014, 2018)⁶ that uses the Extended IRTF Library (E-IRTF; Villaume et al. 2017) and the MIST isochrones (Choi et al. 2016). With alf we can model the full continuum-normalized spectrum of integrated light for stellar ages >1 Gyr and for metallicities in the range of ∼−2.0 to +0.25. The full model has 36 free parameters (see Table 2 in Conroy et al. 2018). The parameter space is explored using a Markov Chain Monte Carlo algorithm (emcee; Foreman-Mackey et al. 2013). In this work we use the priors as described in Conroy et al. (2018) and fix the initial mass function (IMF) to the Kroupa (2001) form.

Theoretical elemental response functions that tabulate the effect on the spectrum of enhancing each individual element modeled in alf were computed with the ATLAS and SYNTHE programs (Kurucz 1970, 1993). For the α-elements relative to Fe considered in our analysis (Mg, Si, and Ca) we correct for the underlying abundance pattern in the empirical stellar library using the [Mg/Fe] values from Milone et al. (2011) and [Ca/Fe] values from Bensby et al. (2014). We assume that Si has the same library abundance pattern as Ca.

We analyze several different data sets in this work (see below). To make the different samples as homogeneous as possible, we fitted over the same spectral range for every spectrum analyzed in this work: 4000 Å < λ < 4400 Å and 4400 Å < λ < 5225 Å.

While we obtain estimates of the light-weighted age as part of the alf models, we do not include age in our analysis. This is because of the uncertain effect of the blue horizontal branch, particularly in the GCs, which can make the inferred ages artificially young. Our analysis of the Milky Way GC system indicates that iron metallicity can still be reliably recovered in the presence of a blue horizontal branch (see Conroy et al. 2018).

Strader et al. (2011) carried out a wide-field kinematic analysis of the M87 GC system using two key data sets: the Keck/LRIS sample and the MMT/Hectospec sample. In Figure 1 we compare the two samples in color–magnitude space. The Keck/LRIS sample (green) was selected to sample the low-luminosity population over the full color range of the GC system, in contrast to previous work that targeted high-luminosity clusters that likely have different properties from the bulk of the GCs (see discussion in Villaume et al. 2019). The MMT/Hectospec sample (yellow) was selected from the higher-luminosity population. The MMT/Hectospec objects were also primarily selected at large radii to aid the sky subtraction since Hectospec is a fiber instrument.

Zoom In Zoom Out Reset image size

In Villaume et al. (2019) we applied full-spectrum SPS models to the Keck/LRIS data set of M87 GCs (Strader et al. 2011) to obtain estimates of iron metallicity ([Fe/H]) relative to solar. We refer the readers to the original paper for details on the modeling and validation of the [Fe/H] values for the Keck/LRIS sample. Here, we do the same analysis for the MMT/Hectospec sample. We used the square root of the summed sky spectrum and flux generated by the reduction pipeline as the uncertainty on the individual GCs. The signal-to-noise ratio (S/N) of this data set is in the range of S/N ∼ 1–30 Å⁻¹ with a resolution of 5 Å. The resolution of the data is higher than the native resolution of the models, so we smoothed the data to 200 km s⁻¹ to be consistent with our previous analysis with the Keck/LRIS data.

Before we smoothed, we identified particularly bad sky lines in the spectra at 4040 Å < λ < 4050 Å, 4355 Å < λ < 4365 Å, and 5458 Å < λ < 5470 Å and interpolated over the flux in each spectrum in those wavelength regions. We fitted 156 spectra and rejected 12 spectra from our analysis based on visual inspection of the residuals between the observed spectra and best-fit models. We show successful fits in Figure 2 for a comparatively high S/N spectrum (S/N ∼ 30; brown) and a low-S/N spectrum (S/N ∼ 10; green), with spectral features of particular interest highlighted. The black line and gray band represent the data flux and uncertainty, respectively.

Zoom In Zoom Out Reset image size

For the individual GCs, we focus our analysis on [Fe/H] and summarize our measurements in Table 1. The majority of the GC spectra do not have sufficient S/N to reliably extract more detailed abundance information. In Section 4.3 we describe how we stacked the individual GC spectra and fit the stacks with alf.

Table 1.  Table of Summary Statistics of the [Fe/H] Measurements for the GCs Included in This Work^a

ID	R.A.	Decl.	[Fe/H]	${\sigma }_{[\mathrm{Fe}/{\rm{H}}]}$	Instrument
H47487	187.73553	12.32802	−1.16	0.30	LRIS
H49585	187.67674	12.32961	−0.51	0.41	LRIS
H49328	187.71423	12.32992	−0.78	0.24	LRIS
...
H47487	187.59446	12.02249	−0.80	0.37	Hectospec
H49585	187.52539	12.03362	−0.40	0.23	Hectospec
H49328	187.91104	12.04288	−1.17	0.39	Hectospec

Note.

^aFull [Fe/H] posteriors of all GC measurements can be found on GitHub:github.com/AlexaVillaume/m87-gc-feh-posteriors and archived in Zenodo doi:10.5281/zenodo.3945492.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

We use data from the Mitchell (formerly VIRUS-P) IFU spectrograph at McDonald Observatory (Murphy et al. 2011; spectroscopy obtained via K. Gebhardt 2020, private communication). The S/N of the individual spectra is in the range of ∼20–50 Å⁻¹. We stacked the individual spectra in 10 bins of galactocentric radius by bootstrapping for the median of the individual spectra in a given bin. We used the 50th percentile from the resulting distribution of flux at a given wavelength as the stacked spectrum and used the average of the 16th and 84th percentiles as the uncertainties on the stacks with the S/N of the stacks ranging from ∼40 to 200, with the outermost spectrum having the lowest S/N.

In Figure 3 we examine the quality of our fits for spectra in the inner (R_gal ∼ 1.32 kpc, S/N ∼ 200) part of the galaxy and the outer part (R_gal ∼ 19.4 kpc, S/N ∼ 40). We compare the Mitchell spectra (black) with the best-fit model spectrum for the inner region (brown) and the outer region spectrum (green) with selected spectral features highlighted. The gray bands are the flux uncertainty from the data. In Figure 4 we examine the residuals between the best-fit model and input data for all the Mitchell spectra used in this analysis. The residuals are typically small (<5%).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We begin with a model for a single population of objects as a way to demonstrate some of the key reasons for using an HBM framework in a simplified setting. Bayesian inference is an application of Bayes's theorem (Bayes & Price 1763),

$$\begin{eqnarray}&&P(A| B)=\displaystyle \frac{P(B| A)P(A)}{P(B)},\ \mathrm{if}\ P(B)\ne 0,\end{eqnarray} \tag{ 1 }$$

which is derived from an axiom of conditional probability. Bayes's theorem is just a way to compute conditional probabilities of events while folding in prior knowledge related to that event. In practice, as a tool for statistical inference, Bayes's theorem is often written in terms of parameters, θ, and data, x, and the denominator, also known as the Bayesian evidence, is often dropped to yield the unnormalized posterior density, $p(\theta | x)\propto p(x| \theta )p(\theta )$ (Gelman et al. 2013).

The first term on the right-hand side of the proportionality is the likelihood function, and the second term is the prior distribution. The likelihood function describes the connection between the available data and the parameters of interest.

In this work, the data we have are $x=[{r}_{\perp ,\mathrm{obs},n},{Z}_{\mathrm{obs},n},{\sigma }_{{\rm{Z}},n}]$ for each n GCs, and our ultimate parameters of interest are the slope m and intercept b (highlighted with a red node and together in the same node to indicate their covariance) of the metallicity gradient. However, we construct our model based on the idea that the data we have correspond to true versions of the parameters that are obfuscated by systematic and random uncertainty. This introduces latent, i.e., unobserved, parameters to our model. That is, the observed metallicity of a GC, ${Z}_{\mathrm{obs},{\rm{n}}}$, is a noisy realization of that GC's true metallicity, ${Z}_{\mathrm{true},n}$, and its observed projected distance, ${r}_{\perp ,\mathrm{obs},n}$, is a realization of the true 3D distance, $| {{\boldsymbol{r}}}_{n}|$.

On the left-hand side of Figure 5, we show part of the graphical representation of our probabilistic model ($| {{\boldsymbol{r}}}_{n}|$ will be discussed in more detail later). This shows the relationship between the observations (filled nodes) and the parameters relevant to the inference we want to make (open nodes). Within the rectangle (known as the "plate"), we show the data and parameters for the individual GCs in the sample. The parameters outside the plate are the parameters for the whole population. We will now distinguish these population parameters (the hyperparameters, α) from the parameters for the individual GCs (θ_n). With the introduction of the α parameters, the joint distribution we seek to constrain is $p(\alpha ,{\theta }_{n}| {x}_{n})$, to which we can apply Bayes's theorem:

$$\begin{eqnarray}&&p(\alpha ,{\theta }_{n}| {x}_{n})\propto p({x}_{n}| {\theta }_{n},\alpha )p({\theta }_{n},\alpha ).\end{eqnarray} \tag{ 2 }$$

Zoom In Zoom Out Reset image size

In Figure 5 the arrows represent the conditional dependency among the different parameters and makes clear the hierarchical nature of our model. The key point is that the data, x_n, are only conditionally dependent on the parameters θ_n and are therefore conditionally independent from the hyperparameters, α. That, as well as being able to factor p(θ_n, α) to $p({\theta }_{n}| \alpha )p(\alpha )$, gives

$$\begin{eqnarray}&&p(\alpha ,{\theta }_{n}| {x}_{n})\propto p({x}_{n}| {\theta }_{n})p({\theta }_{n}| \alpha )p(\alpha ).\end{eqnarray} \tag{ 3 }$$

The key difference between standard Bayesian models and HBMs is that we constrain the population parameters by conditioning on the observations of the many individual GCs rather than fixing them and using them as priors.

The gradient parameters are inferred through modeling the "true" metallicity values of the individual GCs, ${Z}_{\mathrm{true},n}$. The ${Z}_{\mathrm{true},n}$ values are set deterministically by the linear relation p(${Z}_{\mathrm{true},n}| m,b,{{\boldsymbol{r}}}_{n})=m\times | {{\boldsymbol{r}}}_{n}| +b$. We condition ${Z}_{\mathrm{true},n}$ on ${Z}_{\mathrm{obs},n}$ by modeling the observations as drawn from normals centered on the true values and with a standard deviation that encompasses our uncertainty on the metallicity gradient. This uncertainty is the quadrature sum of the uncertainties on the individual [Fe/H] measurements, ${\sigma }_{Z,n}$, and unobserved uncertainty for the intrinsic scatter in the radial metallicity gradient, S, such that $\sigma =\sqrt{{S}^{2}+{\sigma }_{Z,n}^{2}}$.

With the ${{\boldsymbol{r}}}_{n}$ dependence for ${Z}_{\mathrm{true},n}$ we introduce a key advantage when using a Bayesian framework. Even though we do not have the line-of-sight distances, ${r}_{\parallel ,n}$, we can make inferences on the true distances for each GC while only making weak assumptions about the population. Specifically, we model the angular distribution of GCs as isotropic and assume that the GCs are normally distributed in radius by some scale length, R, in all three coordinates xyz and marginalize over the angle, ϕ_n, between x_n and y_n to get

$$\begin{eqnarray}\begin{array}{rcl}p({{\boldsymbol{r}}}_{n}| R) & = & \displaystyle \frac{{r}_{\perp ,n}}{{R}^{2}}\exp \left[\displaystyle \frac{-{r}_{\perp ,n}^{2}}{2{R}^{2}}\right]\\ & & \times \,\displaystyle \frac{1}{\sqrt{2\pi {R}^{2}}}\exp \left[\displaystyle \frac{-{r}_{\parallel ,n}^{2}}{2{R}^{2}}\right],\end{array}\end{eqnarray} \tag{ 4 }$$

which we fully derive in Appendix A. As such, we model the projected distances as drawn from a Rayleigh distribution (the first term on the right-hand side of the above equation) and the line-of-sight distances as drawn from a normal distribution (the second term). This structure is graphically represented on the right-hand side of Figure 5.

In reality, a power-law distribution better describes the projected radial distribution of a typical GC system. In the left panel of Figure 6, we compare the expected quantiles from a Rayleigh distribution versus the quantiles of the projected distances for our M87 sample (open circles). This figure demonstrates the extent the mock spatial distribution deviates from the assumption of our model. If the projected distances were drawn from a Rayleigh distribution, the theoretical quantiles versus the data quantiles would be a straight line. The Rayleigh distribution is not a significant deviation from the distribution of the observed projected distances in our sample.

Zoom In Zoom Out Reset image size

Figure 5 displays the joint probability distribution of all our parameters and data, $p({Z}_{\mathrm{obs},n},{r}_{\perp ,\mathrm{obs},n},{Z}_{\mathrm{true},n},{{\boldsymbol{r}}}_{n},R,S,m,b,{\sigma }_{n})$. Because the arrows indicate the conditional dependence among the parameters and data, we use this to factorize the joint probability distribution of all our parameters into conditionally independent probability distributions to obtain

$$\begin{eqnarray}&&\begin{array}{l}p({Z}_{\mathrm{obs},n},{r}_{\perp ,\mathrm{obs},n},{Z}_{\mathrm{true},n},{{\boldsymbol{r}}}_{n},R,S,m,b,{\sigma }_{n})\\ \propto \ \displaystyle \prod _{n=1}^{N}p({Z}_{\mathrm{obs},{\rm{n}}}| {Z}_{\mathrm{true},{\rm{n}}},S,{\sigma }_{n})p({r}_{\perp ,{obs},n}| {{\boldsymbol{r}}}_{n})\\ \quad \times \ \displaystyle \prod _{n=1}^{N}p({Z}_{\mathrm{true},n}| m,b,{{\boldsymbol{r}}}_{n})p({{\boldsymbol{r}}}_{n}| R)p({\sigma }_{n})\\ \quad \times \ p(S)p(R)p(m,b).\end{array}\end{eqnarray} \tag{ 5 }$$

To test the efficacy of this model, we generated mock data from where the coordinates are drawn from a power law that goes as −2.5, and each data point is randomly assigned 10%−55% uncertainty. In the right panel of Figure 6 we compare how well we recover the true slope when using projected distance as a proxy for true distance and a weighted least-squares fit to get the gradients (open circles) to how well we recover the slope when we marginalize over the scale length (closed circles). Over the range of slope values, the recoverability improves with the statistical deprojection. The difference in results between the two methods is starkest when the gradient is significant, while there is no difference in the recoverability when the gradient is consistent with being flat (m = 0).

3.2.1. Effect of Making Cuts on the Population

In the previous section, we demonstrated the efficacy of a Bayesian linear regression approach relative to weighted least-squares to accurately recover the gradient parameters of a set of data points drawn from a particular line. A fundamental assumption in the method presented is that the data points come from the same population. As previously described, however, GC systems are generally composed of subpopulations, and knowing how to separate the individual GCs of a system into the correct subpopulations is one of the most difficult steps toward characterizing GC systems.

In the top panels of Figures 7–9 we show three versions of mock data, all generated from power-law distributions and two underlying gradients. In all versions we use ${b}_{\mathrm{true},0}=-0.4$ and ${b}_{\mathrm{true},1}=-1.0$ and a variety of slope parameters: ${m}_{\mathrm{true},0}\,={m}_{\mathrm{true},1}=-0.05$ (Figure 7), ${m}_{\mathrm{true},0}={m}_{\mathrm{true},1}=-0.01$ (Figure 8), and ${m}_{\mathrm{true},0}\,=-0.01$ and ${m}_{\mathrm{true},1}=-0.015$ (Figure 9). For each "system" we generated five realizations of mock data.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For each set of mock data, we made constant cuts based on the MDFs to separate the populations, mimicking what one might do if one did not have a priori information on the different subpopulations. We fit each realization of the subsequent subpopulations with our Bayesian linear regression model presented in Section 3.1, which we note is already an improvement over previously used methods, as demonstrated in Section 3.1.

In the middle panels of Figures 7–9 we show how effective this method is by comparing the true slope values (black line) to the median of the posteriors for each realization of the mock data (brown histograms). Even with the improvements to the linear regression outlined in Section 3.1, the inferred slope values are not accurate, with the inferred slopes typically being flatter than the true slopes. We therefore need to generalize our single-population model to account for the covariance between the gradient parameters and the subpopulation membership assignments to more accurately estimate both.

3.2.2. A Mixture Model

Hogg et al. (2010) discussed mixture models in the context of linear regression for the purposes of outlier rejection. Separating individual GCs into subpopulations is an equivalent problem. We model the system such that a given GC has C number of subpopulations it could be assigned to through an identifier parameter q_n. Like Hogg et al. (2010), we directly marginalize over the class membership of each GC by introducing a new parameter, the prior on q_n, P_c ∈ [0, 1] such that ${\sum }_{c=1}^{C}{P}_{c}=1$. The P_c parameters are the mixture weights for each subpopulation and allow us to marginalize out the subpopulation identifiers.

In principle, we can fit for any number of subpopulations. In practice, however, throughout this work we specialize to the C = 2 case such that we model the mock and observed data as a bimodal distribution. We set a lower limit on P_c, P_min = 0.3. Then,

$$\begin{eqnarray}&&{P}_{c}=\left\{\begin{array}{l}{P}_{0}\sim \mathrm{Uniform}({P}_{\min },1-{P}_{\min })\\ {P}_{1}\sim 1-{P}_{0}\end{array}\right..\end{eqnarray} \tag{ 6 }$$

The structure otherwise remains the same as our single-population model. We are able to transition our population parameters from the single-population model to be a part of the mixture model because the parameters will exist for each mixture component (i.e., GC subpopulation). So we make a small adjustment to our notation: m_c, b_c, log R_c, and log S_c, where the subscript refers to a given subpopulation. The joint probability distribution is then given by

$$\begin{eqnarray}&&\begin{array}{l}p({\alpha }_{C},{\theta }_{n}| {x}_{n})\propto \left(\displaystyle \prod _{c=1}^{C}p({R}_{c},{m}_{c}{b}_{c},{S}_{c})\right)\times \displaystyle \prod _{n=1}^{N}\\ \ \times \,\left(\displaystyle \sum _{c=1}^{C}{P}_{c}\times p({Z}_{\mathrm{obs},n},{r}_{\perp ,\mathrm{obs},n},{Z}_{\mathrm{true},n},{{\boldsymbol{r}}}_{n},{R}_{c},{S}_{c},{m}_{c},{b}_{c},{\sigma }_{n})\right).\end{array}\end{eqnarray} \tag{ 7 }$$

The third term in this equation is what we factorized in Section 3.1 for the single-population model. We show the graphical representation of the final hierarchical mixture model in Figure 10.

Zoom In Zoom Out Reset image size

We specify our model with the probabilistic programming package PyMC3 (Salvatier et al. 2016). PyMC3 uses the Hamiltonian Monte Carlo (HMC) family of samplers. For this work in particular, we use the No-U-Turn Sampler (NUTS; Hoffman & Gelman 2011). HMC samplers are more efficient than the commonly used ensemble samplers because they do not rely on the current state to propose the next state (for an introduction to HMC see Betancourt 2017), and so they are the more appropriate choice for high-dimensional problems.

The sampling efficiency of the HMC algorithm is highly sensitive to several tuning parameters. For this work, the most important tuning parameter is the mass matrix because our model parameters are highly covariant. If the mass matrix is not well matched to the covariance of the posterior, both the step size will need to be decreased and the number of steps will need to be increased, making it difficult to achieve convergence.

PyMC3 does not have a built-in way to optimize the mass matrix. We use the exoplanet⁷ extensions to PyMC3 to fit for a dense mass matrix during burn-in. We find values to initialize the sampler via several steps: first, we fit a 1D mixture of Guassians on the metallicities while taking into account the uncertainties to get an initial guess of the class membership for each observation, and second, we fit a linear model to the project metallicity gradients for each subpopulation to find initial guesses for the intercepts and slopes. With this approach, we obtain a converged model based on the Gelman-Rubin statistic for each parameter, $\hat{R}$ (where $\hat{R}\gt 1$ indicates that the chains have not converged).

In the bottom panels of Figures 7–9 we show the inferred slope posteriors (black histograms) to demonstrate the efficacy of this method. In all cases the recovery is better when using the full model, with the biggest improvement made in the case where the two subpopulations are most well mixed in the MDF (Figure 7).

The full model accurately recovers the different slopes in Figure 9 within 1σ uncertainty but cannot distinguish the gradients as different at a statistically significant level. This is still an improvement over existing methods, but, in the context of GC subpopulations, the ability to discern any gradient differences is essential for understanding how potentially similar the assembly histories of the different subpopulations are (see Section 5 for more discussion on this). Improving the precision of the subpopulation parameters will be the subject of future work.

In this work, we have two spectroscopic data sets for the GC system: the Keck/LRIS sample and the MMT/Hectospec sample. The former covers a radial range of ∼7–27 kpc, while the latter spans ∼14–142 kpc. Previous kinematical analyses of GCs and planetary nebula show signs of a transition at ∼40–50 kpc, which may be related to a recent accretion event (Romanowsky et al. 2012; Longobardi et al. 2015; Zhang et al. 2015). Photometric surveys have also shown in M87, as well as other massive ETGs, that blue GCs begin to dominate at large radii. To measure the metallicity gradients, we split our sample at 40 kpc. The inner halo sample consists mostly of the Keck/LRIS data, with a small fraction coming from the MMT/Hectospec data. The outer halo sample consists completely of MMT/Hectospec data.

Before discussing the results from fitting the model to the individual [Fe/H] measurements, we first establish our expectations empirically in Figure 11. In the top panel we show two MDFs for the Keck/LRIS sample, one for the inner part of the data set (light-green solid line, r_⊥ ≤ 25 kpc) and one for the outer part (dark-green dashed line, r_⊥ > 25 kpc). There are two distinct peaks in the inner MDF, while the outer MDF has a less significant second, metal-rich peak. In the outer bin, the metal-poor peak shifts noticeably from the inner metal-poor peak.

Zoom In Zoom Out Reset image size

In the bottom panel we do the same demonstration for the outer halo. Bimodality is not as clearly seen in the outer halo sample as it is in the inner halo, but there is a distinct negative shift from the main peak from the inner bin to the outer bin. The lack of clear bimodality could be a result of the MMT/Hectospec sample having far fewer red GCs than blue GCs and is consistent with the findings for other BCGs (see, e.g., Harris et al. 2017).

For the modeling, we initialized the MCMC chains in the same manner as the mock data and modeled the data as composed of two subpopulations for both the inner and outer halos. In Figure 12 we show the posteriors of slope values for the metal-poor (blue) and metal-rich populations for the inner halo (top panel) and the outer halo (bottom panel). The 16th and 84th percentiles are marked by the colored bands. In each panel a flat gradient is marked with the black dashed line.

Zoom In Zoom Out Reset image size

The uncertainty on the slope measurements is significantly larger for the inner halo than the outer halo measurements even though the [Fe/H] uncertainty is ∼20% higher for the MMT/Hectospec data. The large uncertainty in the inner halo could be a result of the comparatively nonuniform coverage in r_⊥ for the inner halo sample; we get less information from each individual measurement in the inner halo than the outer halo. For the inner halo, both the metal-rich and metal-poor slopes are consistent with a flat gradient and are statistically consistent with one another. In the outer halo, the metal-poor slope is consistent with a flat gradient, while the metal-rich slope is slightly negative.

In Figure 13 we show the metallicities as a function of the deprojected distances where the data points are colored by subpopulation assignment. The opacity of the individual points is scaled by certainty of that subpopulation assignment, with white indicating that the assignment is highly uncertain. We show the posterior median (solid lines) and the range encompassed by the 16th and 84th percentiles (bands) of the gradient distributions. Even though the gradient parameters are more uncertain in the inner halo, the subpopulation membership assignments are more certain than in the outer halo population because there are fewer metal-rich GCs and the metallicity separation between the subpopulations is smaller.

Zoom In Zoom Out Reset image size

In Figure 14 we compare the MDFs for the metal-poor GCs (blue) and metal-rich GCs (red) for the inner halo GCs (left) and outer halo GCs (right). The solid lines show the result of using the posterior median of the subpopulation assignments of the individual GCs. We represent how the uncertainty in the subpopulation assignment affects the MDF by plotting the result of selecting class labels using a random number generated by the probability of the cluster–subpopulation pair for 10 random samples from the posterior (filled histograms).

Zoom In Zoom Out Reset image size

To check the results of our model, we compare Figure 11 with Figure 14. Figure 11 indicates that the two subpopulations should be of about equal size for the inner halo sample and that the metal-poor GCs would be a larger population in the outer halo sample. Even though the subpopulation membership assignments are much less certain for the outer halo sample, we see that the model assigns significantly fewer metal-rich GCs in the outer halo. This picture is overall consistent with our broad understanding that with increasing galactocentric distance there will be more metal-poor GCs relative to metal-rich ones.

In Table 3 we summarize the gradient parameters and the characteristics of the MDFs for the subpopulations.

In the Milky Way GC system, the metallicity subpopulations are associated with different spatial and kinematical components of the Galaxy itself (Zinn 1985). In Figure 15 we show how well this pattern holds for M87 by examining the chemodynamics of the subpopulations in [Fe/H]–radial velocity space for the inner halo (top) and outer halo (bottom). For the inner halo, our modeled metallicity subpopulations correspond to differences in the radial velocity distributions. The metal-rich subpopulation has significantly less radial velocity dispersion than the metal-poor subpopulation.

Zoom In Zoom Out Reset image size

Unlike the inner halo, there does not seem to be a correspondence between metallicity subpopulation and differences in the kinematic properties of the subpopulations for the outer halo. The subpopulations in the outer halo have a radial velocity dispersion similar to the inner halo metal-rich subpopulation.

In Figure 16 we show the stellar population radial profiles for M87 from our fits to the Murphy et al. (2011) sample (black circles) for [Fe/H] (top left) and a variety of α-elements. The M87 starlight shows a declining [Fe/H] profile and slightly positive profiles for [Mg/Fe] (top right), [Si/Fe] (bottom left), and [Ca/Fe] (bottom right). The [Fe/H] gradients are consistent with previous works that have studied stellar population gradients in massive ETGs (e.g., Greene et al. 2015; van Dokkum et al. 2017b; Gu et al. 2018b). These previous studies typically found a flatter [Mg/Fe] than what we present here, but it is not a substantial difference.

Zoom In Zoom Out Reset image size

In the top left panel we also show [Fe/H] estimates from deep broadband photometry of resolved stars in M87 (contours; Bird et al. 2010). This comparison shows a metal-poor field star population that is not probed by the stellar population models. The comparison of the resolved stars to the integrated light demonstrates the limitations inherent in studies using integrated light and the need for the inclusion of the GC population in the analysis.

To obtain abundance information for the GCs, we have to stack the individual spectra since the majority of the Keck/LRIS and MMT/Hectospec spectra have too low S/N to reliably extract abundance information. Stacking the GC spectra is made difficult by the need to separate the sample by subpopulation, as it is expected that the different subpopulations will have different origins and thus different abundance patterns. We demonstrated in Section 3 the importance of using an HBM framework for making accurate determinations of subpopulation membership of the individual GCs, which help make more physically appropriate stacks. Additionally, we can take advantage of having made probabilistic determinations of subpopulation membership for the individual GCs. In Figure 14 we demonstrated how uncertainty in the subpopulation memberships propagated to the MDF of the GC system. In the same manner, we can propagate that uncertainty to our stacks and abundance information.

In the same manner we used for the Mitchell data, we made four inner halo stacks, binning by metal-rich and metal-poor and then further separating the GCs at a radius at 16 kpc, and two stacks for the outer halo only binning by metal-rich and metal-poor. We made 10 iterations for each subpopulation, determined by different draws from the posterior for different subpopulation assignments for each GC (same draws that are shown in Figure 14). The stacked GC spectra have a typical S/N of ∼100–150/Å.

Each version of each subpopulation stack was fitted using alf in the same manner as the Mitchell data. For each parameter of interest, we computed the 16th, 50th, and 84th percentiles of the posteriors for each fit. In Figure 16 we show the results for the metal-rich stacks (red circles) and metal-poor stacks (blue squares). The inner halo stacks are open symbols, and the outer halo stacks are filled symbols.

While we note that Figure 16 cannot be directly compared to Figure 13 because we have moved from deprojected to projected distances, broadly the [Fe/H] gradient measured from the metal-rich inner halo stacks is consistent with the flat gradient measured from the individual [Fe/H] measurements. For the inner halo metal-poor stacks we find a slightly negative gradient from the stack measurements that differs from the flat gradient shown in Figure 13. The likely cause of this difference is the nonuniform sampling of the inner halo GCs in galactocentric radius. For the individual [Fe/H] measurements the MMT/Hectospec sample provides the only coverage past ∼27 kpc and the more metal-rich measurements (∼ − 1.0) seem to be enough to keep the gradient nearly flat. However, for the stacks there are fewer of these comparatively metal-rich GCs than metal-poor ones, so their contribution to the stack is not as important.

The inner halo metal-rich stacks have a similar metallicity to the M87 starlight in the same region, while the metal-poor stacks are less metal-rich by ∼1 dex. For [Mg/Fe] the inner halo metal-rich stacks are less enhanced than the galaxy light, while the metal-poor stacks have similar abundances for the galaxy light. For both outer halo stacks, there is a precipitous drop in Mg enhancement. For [Si/Fe], the inner halo metal-poor stacks are enhanced relative to the galaxy light and metal-rich stacks, which have similar values. There is a noticeable drop in enhancement for the outer halo metal-poor stack, while the outer halo metal-rich stack is consistent with the inner halo measurements. The [Ca/Fe] measurements are consistent between the galaxy light and all the inner halo stacks, and the outer halo stacks are consistent with the inner halo measurements.

In the left panel of Figure 17 we show [Mg/Fe] versus [Fe/H] for the GC stacks and the galaxy data (colors and symbols are the same as in the previous figure). We also show the abundances for the Milky Way stars (purple cloud; from the JINAbase (Abohalima & Frebel 2018); see detailed references in Appendix B), stars in dwarf galaxies around the Milky Way (brown cloud; JINAbase and Bonifacio et al. 2004), and Milky Way GCs fitted from Schiavon et al. (2005) (purple circles; see Villaume et al. 2019, for details on these fits). Also displayed are the abundances for dwarf ellipticals (dEs) in Virgo from two different studies, Şen et al. (2018; open squares) and Sybilska et al. (2018; triangles). For the Sybilska et al. (2018) sample we differentiate between "M87 dEs" (filled) and "Virgo dEs" (open) with a cut at 300 kpc from M87 (distances from Peng et al. 2008). We were motivated by the results from Liu et al. (2016), which showed a transition in [Mg/Fe] in the dwarf elliptical population at this distance.

Zoom In Zoom Out Reset image size

In the right panel of Figure 17 we show the median values for each object of [Si/Fe] and [Ca/Fe], except for Şen et al. (2018), who only measured [Ca/Fe]. The Sybilska et al. (2018) study is not shown in this panel, as they only measured [Mg/Fe]. All symbols are the same as the left panel. The outer halo metal-rich stack is enhanced in [(Si,Ca)/Fe] but much closer to solar in [Mg/Fe].

In Figure 18 we show a similar figure to Figure 16, but now for light elements: radial profiles for [C/Fe] (left) and [N/Fe] (right). For the M87 starlight we see enhanced [C/Fe] and [N/Fe] values and with negative radial gradients for both abundances. The GCs as a whole are less enhanced than the galaxy starlight for [C/Fe] but more enhanced in [N/Fe].

Zoom In Zoom Out Reset image size

Piecing together how the inner halo (<40 kpc) formed is complicated by the fact that it is a mix of in situ and ex situ stellar populations of unknown proportions. Decomposing the whole population into these components from observations cannot be quantitatively done with integrated galaxy starlight alone. With the GC system we have discrete tracers of near-simple stellar populations that overlap with and extend our coverage of the galaxy field star population, providing additional insight into how this region formed.

The red GC populations in massive ETGs have long been thought to have formed along with the original galaxy because they follow the field star density profiles and kinematics (see Strader et al. 2011, for M87 in particular). Until this work, however, a direct metallicity comparison at the same galactocentric radius has not been done. We have established that the metal-rich GCs have an average [Fe/H] ∼ −0.4 and [Mg/Fe] ∼ +0.15 (Figures 16 and 17), similar to the galaxy field star population ([Fe/H] ∼ −0.3 and [Mg/Fe] ∼ +0.40) over the same radial extent, indicating a common origin of the two populations.

In Figure 16 we show the radial gradients for [Fe/H] (top left panel) and various α-elements measured from the integrated light of M87 (black circles). We find that the [Fe/H] gradient for the field star population is flat within the inner 2 kpc and then steepens to a negative gradient. This negative gradient is paralleled by rising gradients in all of the α-elements. This trend is characteristic of the populations seen in other massive galaxies (Greene et al. 2013; Gu et al. 2018a) and can be viewed as consistent with the second phase of the "two-phase" formation framework (Oser et al. 2010) within a preferentially quenched environment (Liu et al. 2016).

It would follow from this scenario that the bulk of the metal-rich GCs came in from mergers. However, the measurements from the GC stacks indicate that the [Mg/Fe] values for the metal-rich GCs decline with radius. The population that then presumably brought in the metal-rich GCs would be diluting the [Mg/Fe] and depressing the gradient, rather than contributing to its rise. On the other hand, the metal-poor GCs are very Mg enhanced.

A negative gradient is expected from the kind of minor mergers that would bring metal-poor GCs into M87. Major mergers, which would bring in metal-rich GCs, can flatten gradients (e.g., Taylor & Kobayashi 2017). Our current measurements show that the metal-poor subpopulation gradient is skewed negative but is consistent with the metal-rich one (and with a flat gradient; Figure 12) within the 1σ uncertainties. It would presumably follow then that the metal-poor GC population was affected by the same processes that flattened the metal-rich GC gradient, i.e., that the metal-poor population was already in place by the time the major mergers began, and thus the low-mass satellites were preferentially accreted early. To the best of our knowledge, this is not something that has been examined in cosmological simulations of galaxy evolution. However, we can use the additional constraint of the metallicity gradient of the galaxy light. If all the low-mass satellite galaxies came in early, we would expect that the stellar metallicity gradient should be similarly flat to the GC subpopulation gradients, and so we conclude that that scenario is not likely.

An alternative explanation is that the metal-poor population is not entirely ex situ. Mandelker et al. (2018) described a scenario in which metal-poor GCs form in situ in massive halos (i.e., viral masses a few times 10¹⁰–10¹¹ M_⊙) at high redshift (z ∼ 6) as a result of fragmentation of the unstable cold gas filaments accreting on the forming galaxy. To the best of our knowledge, at this redshift the circumgalactic gas does not evolve (see Figure 6 in Robert et al. 2019), and so there would be no mass–metallicity relation, and thus the population of GCs that stem from this scenario would have a flat metallicity gradient. An in situ metal-poor population commingling with the ex situ metal-poor GC population could then "dilute" the overall "metal-poor" subpopulation gradient we measure.

A related issue regarding the origins of the inner halo metal-poor GC population is that cosmological simulations of galaxy evolution do not predict such a low-metallicity population. Cosmological simulations predict that the Milky Way–mass (mass ratio ∼1:5) galaxies are the primary building blocks of the stellar halos of massive galaxies (Oser et al. 2012; Pillepich et al. 2018), while the typical mass implied by the median metallicity of the metal-poor GC population is much lower, M_* ∼ 10⁸ (Kirby et al. 2013). This makes the metal-poor GCs too metal-poor to fit this framework, even though in Villaume et al. (2019) we established that the metal-poor GCs in the inner halo of M87 are ∼0.4 dex more metal-rich than previously thought. The mass ratios of the mergers suggested by this metal-poor population are not inconsistent with analytic models of merger-driven BCG expansion from high redshift to the present day (see Equation (4) in Naab et al. 2009).

5.1.1. Metal Production in the Field and GC Populations

While observations of GCs and planetary nebulae indicate the presence of an intracluster component that becomes significant beyond ∼300 kpc (Longobardi et al. 2018a, 2018b), our work focuses on the material inside this radius and bound to M87. The presumption is that the stellar population parameters in the outer halo provide cleaner benchmarks by which to judge accretion predictions since a significant in situ population is not expected to complicate the interpretation.

As discussed in the previous section, the metallicity and [Mg/Fe] gradients in the inner halos of M87 and other massive ETGs are evidence that there was preferential accretion of environmentally quenched satellites. This does not appear to be the case for the outer halo, as the outer halo GC stacks have lower [Mg/Fe] abundances than the inner halo stacks. However, when comparing the abundance differences between the inner and outer halo GCs, we need to consider the possibility of mass effects. The outer halo sample consists of more luminous, and therefore more massive, GCs than the inner halo sample (Figure 1). The abundance spreads in the Milky Way GCs have been shown to correlate with luminosity (Figure 16 in Carretta et al. 2010), so it would follow that the more massive outer halo GCs might in some way be impacted by this effect. We need to address whether we expect this effect to be significant.

Carretta et al. (2010) demonstrated that the correlation between luminosity and extent of abundance spreads (specifically in their case, Na–O) in GCs is driven by the extreme of the abundance anticorrelation, not the median values (see their Figure 11). Since integrated light probes the average parameters of the stellar populations, the mass dependency might not be as strong when measuring integrated light. To test this, we found the correlation between mass for the Milky Way GCs and [C/Fe], [N/Fe], and [Mg/Fe] as measured by alf from integrated light. For [Mg/Fe] and [C/Fe] we found a mass dependence of ∼16%. Even accounting for the outer halo GCs being more massive than the most massive GCs, this effect is unlikely to explain the full difference between the inner and outer halo abundances, thereby confirming the lower Mg enhancement in the outer halo compared to the inner halo.

To determine the specifics of the progenitor satellites that built up the outer halo in M87, we find that there appear to be at least two distinct populations in the outer halo, in contrast to previous studies that have treated the outer halo GC populations as monolithic (e.g., Forbes & Remus 2018). First, while the outer halo population has an overall lower metallicity than the inner halo population, the individual GCs display a large range in metallicity. Moreover, even when accounting for the uncertainties in the subpopulation membership assignments when creating the two outer halo stacks, the metallicities of those stacks are different in a statistically significant manner. Second, the metal-poor stack does not display the same unusual α-element abundance pattern as the metal-rich stack or the dEs (see Figure 17).

Even though Mg, Si, and Ca are all α-elements, they have different formation sites. Mg is purely a product of massive stars, while Si and Ca can both also be produced in Type Ia supernovae (Woosley et al. 2002). This opens the possibility of using the abundance patterns of the GC stacks to tag the GC population in M87 to possible progenitors. The α-element abundance pattern displayed by the metal-rich stack is echoed by the dE population in the Virgo Cluster (Figure 17). This makes it tempting to point to the dE population as the progenitors of some portion of the outer halo of M87. However, it is important to acknowledge the difficulty in abundance tagging in this scenario because of the strong radial gradients in dEs (Figure 1; Sybilska et al. 2018). Şen et al. (2018) used the R_e/8 aperture with the nucleus included, while Sybilska et al. (2018) took the luminosity-weighted average of spectra within 1R_e. While noting this caveat, the metallicity gradient for the outer halo metal-rich subpopulation is negative and inconsistent with a flat gradient within the 1σ uncertainty, which is consistent with the bulk of this population being brought in by low-mass satellites.

Furthermore, the progenitor mass implied by the median metallicity of the outer halo metal-rich GC population is M_* ∼ 10⁹ M_⊙ (Gallazzi et al. 2005), which is consistent with the predictions from IllustrisTNG. At >100 kpc, Pillepich et al. (2018) predicted that 90% of the ex situ mass comes from progenitors with stellar masses ≳5 × 10⁹ M_⊙, with the typical progenitor mass being ∼7 × 10¹⁰ M_⊙ (see their Figure 13(b)). However, as in the inner halo, the metal-poor population is not consistent with the predictions from simulations.

Other aspects of the GC system are consistent with our interpretation of the stellar population parameters. In the outer halo, M87 has a V-band luminosity of ∼2.9 × 10¹⁰ L_⊙, very similar to the total luminosity of the Milky Way (Kormendy et al. 2009; Bland-Hawthorn & Gerhard 2016). From Sérsic fits to the photometric sample of M87 GCs from Strader et al. (2011) and correcting for GC luminosity function incompleteness, we estimate that there are ∼1200 metal-rich GCs and ∼4500 metal-poor GCs in this region. The GC specific frequency (S_N) of the outer halo is then S_N ∼ 16; besides M87, the only galaxies in Virgo with such a high value are dwarfs with M_V ∼ −17 to −16 and fainter (${M}_{* }\sim {10}^{8}\mbox{--}{10}^{9}\,{M}_{\odot }$; Peng et al. 2008, Figure 12).

A similar conclusion to our own was reached by Longobardi et al. (2018a) using the outer halo light color M87 to infer low-mass progenitors. Hartke et al. (2018) took that result to indicate a problem with the feedback prescription in IllustrisTNG. However, we also need to consider the possibility that the accreted satellites were different from the surviving population, that is, that there may have been Milky Way–mass galaxies but with S_N more like dwarfs, which no longer exist today, or possibly similar to the ultradiffuse galaxies that have been found to have very high S_N (Peng & Lim 2016). The evidence for a dynamically cold phase-space shell prompted Romanowsky et al. (2012) to suggest that ∼20% of M87's GC population within $50\lt {R}_{\mathrm{gal}}/\mathrm{kpc}\lt 95$ came from an E/S0 progenitor with luminosity ∼0.5L*, which could lead to the flattened metallicity gradient we measure in the outer halo metal-poor GC subpopulation.