Multivariate analysis of the globular clusters in M87

The choice of statistical methods to apply to a set of data depends on its Gaussian or non Gaussian nature. We have thus tested the Gaussianity of our data set by the Shapiro-Wilk test (1965). According to this test the null hypothesis is that the data set is Gaussian. The test statistic W is defined as W=${\Sigma }_{i=1}^n{a}_i{x}_i^2$/${\Sigma }_{i=1}^n(x_i-\overline{x}{)}^2$

where n is the number of observations, $x_i$ are the ordered sample values and $a_i$ are constants generated from means, variances and covariances of the order statistics of a sample of size n from a normal distribution. In our case the data set is multivariate, hence we have used the multivariate extension of the Shapiro-Wilk test. We found that the p value of the test is $1.327\times {10}^{-11}$, which is very small. Thus the null hypothesis is rejected at a 5% level of significance. We conclude that the present data set is non Gaussian in nature.

While Principal Component Analysis (PCA) applies to Gaussian data sets, ICA applies to non Gaussian data sets. ICA is a dimension reduction technique, i.e. it reduces the number of observable parameters p to a number m (m$<<$p) of new parameters, where these new parameters are the linear combinations of p parameters, such that these m parameters are mutually independent. Mathematically speaking let $X_1,X_2,......,X_p$ be p random vectors (here p parameters) and n (here 147) be the number of observations of each $X_i$ (i=1,2,…,p).

Let

where $S=[S_1,S_2,......,S_p{]}^{{}^{\prime }}$ is a random vector of hidden components $S_i$ (i=1,2,…,n) such that $S_i$ are mutually independent and A is a non singular matrix. The goal of ICA is to estimate A and to find S by inverting A i.e.

Under ICA we find the components $[S_1,S_2,......,S_n]$ those are independent by finding an unmixing matrix W in such a way that the covariance between any two nonlinear functions $g_1\left(S_i\right)$ and $g_2\left(S_j\right)$ for $i\ne j$ is zero i.e. the independent components are nonlinearly uncorrelated (for more details see Comon (1994); Chattopadhyay et al. (2013b) and references therein).

At present there is no better method available to automatically determine the optimum number of Independent Components (ICs). In this paper, the number of ICs is determined by the number of Principal Components (PCs) chosen (Albazzaz and Wang, 2004). To reduce the number of components $S_i$ from p to m (m$<<$p), one is to perform PCA (Babu et al. 2009; Chattopadhyay & Chattopadhyay 2007; Fraix Burnet et al. 2010; Chattopadhyay et al. 2010). In this method also $Y_i$’s (i=1,2,…,p) vectors are found, which are linear combinations of $X_i$’s (i=1,2,…,p) such that $Y_i$’s are uncorrelated. The number of (initially p) components is reduced to m by taking those $Y_i$’s (i=1,2,…,m) for which the corresponding eigenvalues ${\lambda }_i\sim 1$.

In the present work we have first performed PCA to find the significant number of IC components to be chosen. PCA applies to Gaussian data, but the present data set is non Gaussian, so we performed ICA. In PCA the maximum variation with significantly high eigenvalue (viz.$\lambda \sim 1$) was found to be almost 78% for three PCA components. Hence we have chosen three IC components for the cluster analysis. The parameters chosen for ICA are the Lick indices $H\beta$, Mg1, Mg2, Mgb, Fe5270, Fe5335, NaD, TiO, the metallicity Fe/H, the photometric parameters magnitude $i_{mag}$ and colours (g-r), (g-i) and the projected galactocentric distance $R$ (in arcsec).

The optimum value of k is found as follows. First one finds groups for k = 1, 2, 3,…etc. Then a distance measure $d_k$ is computed by

which is the distance of $x_k$ vector (values of parameters) from the centroid $c_k$ of the corresponding group. The optimum value of k is that for which $J_k=(d_k^{{}^{\prime }-p/2}-d_{k-1}^{{}^{\prime }-p/2})$, is maximum (Sugar & James 2003). In the present case k = 4.

For cluster analysis (CA) the number of parameters are 3 IC components chosen as previously described. The optimum number of groups is selected objectively by k-means clustering (Mac Queen 1967), together with the method developed by Sugar and James (2003). The optimum number is k = 4. The groups are labelled G1, G2, G3 and G4.

For testing the robustness of the classification we have proceeded in the following manner. In the original scheme the number of parameters is 13. First we have constructed the variance covariance matrix of the parameters and selected the parameters having maximum variances e.g. here $H_{\beta }$, Mgb, Fe5270, Fe5335, NaD, Fe/H, $i_{mag}$, R have variances greater than 0.25 and the remaining ones, Mg1, Mg2, TiO1, g-r, g-i have variances $\sim$0.0. Therefore we have done the analysis with the above 8 parameters and have obtained exactly the same groups with no variation at all. So we can at once say that our classification is robust with respect to the parameters which are responsible for the maximum variation.

Regarding the uncertainties on the various parameters, the method we have used is totally exploratory (i.e. no underlying distribution is assumed). Hence it is not possible to see the effect of error in such type of analysis directly. What we can do at most is to change the values of the parameters within the range permitted by the error bars and redo the analysis. We have tested it for that for few values of the parameters and found the same groups. Alternatively once the optimum classification (clustering) is obtained, one can use a process called discriminant analysis (Johnson & Wichern 1998) to verify the acceptability of the classification of different GCs. In this standard procedure, using the probability density functions in parameter space for the different clusters, one can assign an object (in this case a globular cluster) to be a member of a certain class. If the original classification is robust, then every GC should be classified again as a member of the same class that it was before. If a significant number of objects are not reclassified then that means that the original classification is not stable and hence not trustworthy. Table 1 shows the result of a discriminant analysis, where the columns represent how the GCs of a cluster were assigned by the analysis. The fraction of correct classifications is 0.925, which implies that the classification is indeed robust by computing classification /misclassification probabilities for the GCs.

This algorithm is used to find the most probable mean rotation curve for each group, assuming that they have a solid-body rotation around the centre of the galaxy. This is a reasonable assumption for the two innermost groups, G1 and G4, because the galaxy itself is in solid-body rotation at least to a radius of 225” (Cohen & Ryzhov, 1997). For the two other groups, G2 and G3, we also tried a rotation curve that flattens beyond 225”, but the data were better adjusted by a solid-body rotation curve (also see Fig.7 of Kissler-Patig & Gebhardt 1998).

The rotation amplitudes $v_{rot}$ ( = $\Omega R$) and the position angles (${\psi }_0$) of the axes of rotation (East to North) of the different groups G1 - G4, (found in the CA) are computed by the relation

(Côté et al. 2001; Richtler et al. 2004; Woodley et al. 2007), where $v_r$ is the observed radial velocity, $v_{sys}$ is taken as the recession velocity of the galaxy, which is 1287 km/s.¹¹1Côté et al. (2001) adopted $v_{sys}$=1350 km/s., $R$ is the projected distance of each GC and $\psi$ is its position angle, measured East to North. We have used Levenberg-Marquardt fitting method (Levenberg 1944; Marquardt 1963) to solve for $v_{rot}$ and ${\psi }_0$. They are listed in Table 2 for G1 - G4. Now kinematic data sets for G1 - G4 are small and for such situation Monte Carlo simulations are used to increase the data sets to have a more convincing result. Monte Carlo simulation needs distributional assumption and here no well known bivariate distribution is fitting with the data well. So we have taken several bootstrap samples and have computed the mean values with standard errors of the rotation parameters . The errors are small as seen from Table 2.

The cluster analysis divided the GCs of M87 into four groups. We did not use the kinematical data ($v_{rad}$ and $\psi$) in the statistical analysis. Nor did we use Mgb/Fe. The latter is an indicator of the abundance in light elements. But we did use these parameters for interpreting the results.

The results of the analysis are summarized in Table 2, which lists the mean values of all the parameters with standard errors for the four groups. The mean velocities of rotation and mean position angles have been computed by the Levenberg-Marquardt algorithm.

To show how the groups are separated we have plotted three Independent Components (viz. IC1, IC2 and IC3) for the groups G1 - G4 found in our analysis (Fig.2). That figure shows that the groups are well separated in IC space.

We first investigate how the different groups differ in their Lick index values, and how they compare in this respect with the simple stellar population models of Thomas et al. (2011). As shown on Figs.3, 4 and 5, G3 has the lowest values of Fe5270, NaD,TiO and Mgb, while G4 has the highest values of these Lick indices, and G1 and G2 are in between. These two latter groups differ in that G1 is marginally less chemically evolved than G2 (these two groups also differ in their spatial distribution and in their kinematics). In other words, the order of increasing overall chemical evolution is G3, G1, G2, G4. The Lick index H$\beta$ is an age indicator, higher values meaning younger ages, but this index is also sensitive to the colour of the Horizontal Branch (HB), higher values meaning a bluer HB. A bluer HB in turn can be due to He-enriched stars. So the interpretation of Fig.6, which shows the run of H$\beta$ vs Mgb, is best done with the help of stellar population models, namely those of Thomas et al. (2011). The data are well fit by models of 12Gyr, albeit with a large scatter suggesting a variety of HBs and/or of $\alpha$/Fe. An alternative explanation in terms of stellar populations with younger ages is not possible because it would not be compatible with the observed broad-band colours (see below).

We next investigate the colours of the stellar populations in the GCs of M87. The colour-colour diagram for the different groups is shown in Fig.7. There is a progressive reddening of the populations, from G3 (green) to G1 (blue) followed by G2 (red) and G4 (black). This is predicted by models of synthetic stellar populations as a result of increasing age and/or metallicity. Since the models of Thomas et al. (2011) do not predict the evolution of broad-band colours, we adjusted Yonsei models of simple stellar populations (Chung et al. 2013) to the data. The models for metallicities ranging from [Fe/H] = -2.6 to 0.6, an age of 12 Gyr and $\alpha /Fe$ = 0.3 are shown as solid lines, green and blue for 0 and 70% secondary populations respectively. Changing $\alpha /Fe$ will only shift the model down diagonally on the figure, and changing the initial mass function has little effect on the models (see Chung et al. 2013). The comparison of our data to the Galev models (Kotulla et al. 2009), and to the models of Maraston et al. (2003) gave results (not shown) very close to those of the Yonsei model for 0% secondary populations. For the Maraston models, we converted the CFH-Megcam magnitudes to SDSS using the paper by Betoule et al. (2013).

The colours of the metal-poor GCs (shown on Fig.7) are well fit by the models. On the other hand, the metal-intermediate (G1) and metal-rich (G2) GCs can only be explained by a Helium-enriched secondary stellar population (Chung et al. 2013). No Yonsei model population with a realistic age is able to account for the colours of group G4. We have no explanation for the colors of this group, other than that they are unusual. It is not possible to advocate younger stellar populations instead of He-enriched ones to interpret the colors of the groups, because such populations are bluer in both colours, but more so in (g -i). Further evidence for He enhancement in the metal-intermediate and metal-rich GCs of M87 is the far ultraviolet excess detected in the GCs of M87 by Sohn et al. (2006), an excess which increases with metallicity. Kaviraj et al. (2007) interpret this excess by the possible presence of hot HB stars from super-He-rich stellar populations, even though their predicted ages are unrealistically old.

We next compare the properties of the GCs of M87 with those of other galaxies. The GCs of M87 have values of Fe5270 and H$\beta$ at a given Mgb (see Figs.3 and 6) that are comparable to those of M31, but generally lower values of NaD and TiO at a given Mgb (see Figs.4 and 5), except G4. In the latter two figures, the predictions of the simple stellar population models seem to hold the middle ground between the GCs of M87 and of M31. Concerning NaD, the high values observed in M31 might be partly due to interstellar absorption and thus less reliable. We are indebted to the referee for this remark. To better understand the significance of these lower values, we also plotted the same data for GCs of two other elliptical galaxies in Virgo, NGC 4472 and 4636 on Figs.4 and 5. The GCs of these two other ellipticals show the same behaviour as those of M87, suggesting that GCs in elliptical galaxies are different from those of spirals in this respect, except at very high metallicities. We speculate that Na is underabundant in G1-G3 and that the Na - O anti-correlation present in Galactic GCs, interpreted as a result of multiple stellar populations, is only present in G4.

The abundance of light elements tends to decrease with increasing metallicity in the GCs of our Galaxy (Fraix-Burnet et al. 2009), M31 (Colucci et al. 2012a) and the LMC (Colucci et al. 2012b), as expected from arguments of nucleosynthesis. High values of $\alpha$/Fe result from star formation that occurs before type 1a supernovae significantly increase the abundance of Fe in the interstellar medium, and are expected mainly in metal-poor environments. In M87, we find that the opposite is true. We use Mgb/Fe = 2Mgb/(Fe5270 + Fe5335) as an indicator of $\alpha$/Fe and calibrate this indicator with the stellar population models of Thomas et al. (2011). As shown on Fig.9, the light-element abundance increases with metallicity. For G1 and G2 it increases from 0 to 0.3, and for G4 it is about 0.4. For G3, which is very metal-poor, the disagreement between $\alpha$/Fe and its indicator Mgb/Fe suggests that the light-element ratio estimates, which are more uncertain at low metallicity, are just not reliable for this group, which is ignored. The contrast between the GCs of M87 and those of our Galaxy and M31 in this respect are further evidence that GCs of elliptical and spiral galaxies have different chemical histories.

We have adjusted a mean rotation curve to each of the four groups (see Sect.3.4). The resulting parameters are listed in Table 2, and the mean rotation curves of the four groups along with radial velocities of the GCs are shown in Fig.10. The only groups with clear evidence of rotation are G2 and G4. The outer (radius $>$ 300 arcsec) clusters of the group G2 rotate faster than the inner ones. Schubert et al.(2010) have found that the metal poor GCs in the outer region of NGC 1399 (the non-rotating giant elliptical galaxy in the centre of the Fornax cluster) is the only GC population that shows significant rotation. The difference with M87 is that G2 is not the metal poorest group. There is marginal evidence for rotation in G1, about an axis which is orthogonal to that of G2. There is also evidence for rotation in G4, whose spatial distribution is flattened in the plane of rotation (along position angles of maximum or minimum radial velocity; see Fig.10). As G4 is centrally concentrated, it is worth mentioning the recent study of Emsellem et al. (2014), who find a counterrotating core in the central 2 arcsec of M87, in position angle 20 degrees. In a detailed kinematic study of the GCs of M87, Strader et al. (2011) have found subpopulations with distinct kinematical behaviours, and concluded that M87 is still in active assembly. They also warned that the old kinematical data (which include the present ones) predict a higher rotation than their new data. Our own kinematical results should thus be viewed with caution.

The first stage of formation of a galaxy is a dissipative collapse. During that phase, chemical enrichment is expected to be accompanied by spin up of the metal-richer material which ends up in circular orbits, whereas the metal-poor material is in eccentric orbits (Eggen et al. 1962). The formation of G3 and G4 can be associated to that phase : G3 in the outer parts from metal-poor material, G4 in the inner part from metal-rich material. The lack of rotation in G3, the flattened spatial distribution, strong rotation and high average $\alpha$/Fe in G4, are properties consistent with this assumption. Montes et al.(2014a) have shown that the central part of M87 is super solar and presumably formed first during the monolithic collapse of the galaxy.

Montes et al.(2014b) have found that the innermost GCs of M87 (within 30 arcsec) are metal poorer than the stars in the region. We show that there are actually two populations of GCs (G3 and G4) near the centre, one of solar metallicity like the stars of the central region (Montes et al. 2014a).

In the second stage of formation, the galaxy accretes low mass satellites in a dissipationless fashion (Bekki et al. 2005). The GCs of G1 and G2 may have formed in this phase. The difference between G1 and G2 could be one of time, G2 forming first and retaining the angular momentum of the gas from which it formed and G1 forming later from cooling flows close to the centre, the gas having lost its angular momentum to collision with gas clouds (viz. Table2, $v_{rot}$). This would explain why the GCs of G1, near the centre, have little rotation, while those of G2, in the outer regions, rotate significantly. In this scenario, one thus expects the velocity of rotation to increase outward, which seems to be the case for the GCs of G2 : $v_{rot}\sim$ 374.93 km$s^{-1}$ for $R\ge$ 300, and $v_{rot}\sim$ 179.92 km$s^{-1}$) for $R<$ 300.