Multivariate analysis of the globular clusters in M87

Abstract: An objective classification of 147 globular clusters in the inner region of the giant elliptical galaxy M87 is carried out with the help of two methods of multivariate analysis. First independent component analysis is used to determine a set of independent variables that are linear combinations of various observed parameters (mostly Lick indices) of the globular clusters. Next Kmeans cluster analysis is applied on the independent components, to find the optimum number of homogeneous groups having an underlying structure. The properties of the four groups of globular clusters thus uncovered are used to explain the formation mechanism of the host galaxy. It is suggested that M87 formed in two successive phases. First a monolithic collapse, which gave rise to an inner group of metalrich clusters with little systematic rotation and an outer group of metalpoor clusters in eccentric orbits. In a second phase, the galaxy accreted lowmass satellites in a dissipationless fashion, from the gas of which the two other groups of globular clusters formed. Evidence is given for a blue stellar population in the more metal rich clusters, which we interpret by Helium enrichment. Finally, it is found that the clusters of M87 differ in some of their chemical properties (NaD, TiO1, light element abundances) from globular clusters in our Galaxy and M31.
Keywords: elliptical galaxies, globular star clusters, data analysis, statistical
1 Introduction
The details of galaxy formation and evolution remain
one of the great unsolved problem in modern astrophysics. Giant
elliptical galaxies at the centre of galaxy clusters are of strong
interest as they provide case studies for galaxy formation theory.
On the other hand globular clusters (GCs) are excellent tracers of
galaxy halos. Indeed some 50% of the stars in our Galaxy may have
originated in GCs (Martell & Grebel 2010). GCs are known to form
in vigorous star formation events (Brodie & Strader 2006 and
references therein) and were an important mode of star formation
in the early universe (Muratov & Gredin 2010). It is possible
that all stars may have originally formed in star clusters which
are only later dispersed into the smooth stellar halos we see
today (Lada & Lada 2003). For making progress in galaxy formation
theory, a holistic approach is required, which involves
compositional as well as kinematical properties of GCs in giant
elliptical galaxies. A landmark study (Emsellem et al. 2007)
suggested that a fundamental dichotomy exists between dynamically
’slow’ and ’fast’ rotating elliptical galaxies. These observations
prompted new interpretations and theoretical models to explain
this difference (Krajnovic et al. 2008; Jesseit et al. 2009).
Proctor et al. (2009) suggested that the rotator class may
change when the kinematics are probed beyond the inner regions.
According to various studies, classical formation of
galaxies have been proposed along five major lines: (i) the
monolithic collapse model, (ii) the major merger model, (iii) the
multiphase dissipational collapse model, (iv) the dissipationless
merger model and (v) accretion and in situ hierarchical merging.
But no model is globally
acceptable over others.
In this context, one is tempted to make a detailed
study of an archetype of each category of elliptical galaxy. In
previous papers (Chattopadhyay et al. 2009; Chattopadhyay et al.
2013b) we studied the GCs of NGC5128, which is a slowly rotating
elliptical galaxy. In the present paper we study the GCs of M87,
which is a strongly rotating elliptical galaxy, located at the
centre of the Virgo cluster. It has a large population of GCs with
wellknown kinematic (Huchra & Brodie 1987; Cohen & Ryzhov 1997;
Hanes et al 2001, Côté et al. 2001; Strader et al. 2011),
photometric (Strom et al. 1981; Strader et al. 2011) and chemical
properties (Mould et al. 1987, 1990; Cohen et al. 1998). To this
end, we carry out a multivariate analysis of the Lick indices,
metallicities and radial distances of the GCs. First we use
independent component analysis (hereafter ICA), suitable for
nonGaussian data sets to search for the independent components
(hereafter IC). The number of such components, which are linear
combinations of parameters, is equal to the number of parameters
(p, say). Since only a few (say, m) of the IC components may
explain a larger percentage of variation in the data, one can take
only those m components instead of all p components. Then the GCs
are classified on the basis of those m independent components
using another exploratory data analytic method, namely Kmeans
cluster analysis (Chattopadhyay et al. 2009, 2010, 2012, 2013a,b)
to find the homogeneous groups. In the end, the properties of the
homogeneous groups of GCs allow us to propose a possible scenario
for the formation of
the GCs and their host galaxy.
The data sets used are presented in section 2. Section 3 describes different methods used in the study. The results and interpretations have been included in Section 4.
2 Data set
The Lick indices (, Mg1, Mg2, Mgb, Fe5270,
Fe5335, NaD and TiO1) and metallicities (Fe/H) used in the
present analysis were taken from Cohen et al. (1998). From the
original sample of 150 GCs we removed the GCs with identification
numbers 5024 and 5026, because they are redundant with 978 (Hanes
et al. 2001), as well as 321 which is a star (Strader et al.
2011), thus reducing the sample to 147 GCs. In the latter paper
the identification number given by Cohen et al. (1998) is preceded
by ”S”. The velocities , halflight radii ,
magnitudes () and and colours were
taken from Strader et al. (2011). The radial distances and
position angles were derived from the coordinates listed in
Strader et al. (2011), assuming that the centre of of the system
of GCs is also that of the galaxy, which is at = 12h
30m 49.42s
and = 12d 23’ 28.044” (Lambert & Gontier 2009).
It is to be noted that Our sample size is small compared
to the recently observed samples by various authors (e.g. Strader
et al. 2011). We have thus plotted the histograms of i 
magnitudes and radial distances of our sample and have compared
with the sample used by Strader et al. (2011) (Fig.1). Now Strader
et al.(2011) have found that the bright and faint GCs have the
same velocity dispersions (sigma = 340 and 339 km /s). In their
section 6.3.4, they find that there is no strong trend with
magnitude for the red GCs. For the blue GCs there is perhaps an
avoidance of the systematic velocity, although marginal. The
conclusion is that the bright and faint compact blue GCs belong to
the same population. So we can safely say that, if the kinematics
are an indicator of the origin of GCs, then our sample
should be enough to determine the origin of the inner GCs.
For comparison purposes, we used the Lick indices of 313 old clusters in M31 (Schiavon et al. 2012), those of 47 GCs in 4472 (Cohen et al. 2003) and 33 GCs in NGC 4636 (Park et al. 2012).
3 Statistical methods
3.1 ShapiroWilk test
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 ShapiroWilk test
(1965). According to this test the null hypothesis is that the
data set is Gaussian. The test statistic W is defined as
W=/
where n is the number of observations, are the
ordered sample values and 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 ShapiroWilk test. We found that the p value of the test is
, 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.
3.2 Independent component analysis
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 (mp) 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 be p random
vectors (here p parameters) and n (here 147) be the number of
observations of each
(i=1,2,…,p).
Let
(1) 
where is a random vector of hidden components (i=1,2,…,n) such that 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.
(2) 
Under ICA we find the components those are independent by finding an unmixing matrix W in
such a way that the covariance between any two nonlinear functions
and for 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 from p to m
(mp), 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 ’s (i=1,2,…,p) vectors are
found, which are linear combinations of ’s (i=1,2,…,p) such
that ’s are uncorrelated. The number of (initially p)
components is reduced to m by taking those ’s (i=1,2,…,m)
for which the corresponding eigenvalues .
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.) 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 , Mg1, Mg2, Mgb, Fe5270, Fe5335, NaD, TiO, the metallicity Fe/H, the photometric parameters magnitude and colours (gr), (gi) and the projected galactocentric distance (in arcsec).
3.3 Kmeans cluster analysis
Kmeans cluster analysis (CA) is a multivariate technique for finding homogeneous groups in a data set giving information of the underlying structure in this data set. In this method one finds k groups, provided each group contains an object and each object belongs to exactly one group. Hence the maximum and minimum number of possible groups are k = n and k = 1 respectively. The algorithm is as follows :

(i) All the objects are divided into k (given) groups in a random manner.

(ii) Any group is selected and any particular object of that group is taken first. Then the parametric distance (here the parameters are 3 IC components) from the chosen object is computed for the remaining objects. If the distance between objects in the group is greater than that for objects in other groups, they are interchanged.

(iii) Process (ii) is applied for all objects and for all groups.

(iv) Steps (i) and (iii) are continued until there is no further change.
The k groups thus found are coherent in nature.
The optimum value of k is found as follows. First one finds groups for k = 1, 2, 3,…etc. Then a distance measure is computed by
(3) 
which is the distance of vector (values of
parameters) from the centroid of the corresponding group.
The optimum value of k is that for which , 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 kmeans 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 , Mgb, Fe5270,
Fe5335, NaD, Fe/H, , R have variances greater than 0.25
and the remaining ones, Mg1, Mg2, TiO1, gr, gi have variances
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.
3.4 LevenbergMarquardt Algorithm
This algorithm is used to find the most probable mean
rotation curve for each group, assuming that they have a
solidbody 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 solidbody 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 solidbody rotation curve (also see Fig.7 of KisslerPatig & Gebhardt 1998).
The rotation amplitudes ( = ) and the position angles () of the axes of rotation (East to North) of the different groups G1  G4, (found in the CA) are computed by the relation
(4) 
(Côté et al. 2001; Richtler et al. 2004; Woodley et al. 2007), where is the observed radial velocity, is taken as the recession velocity of the galaxy, which is 1287 km/s.^{1}^{1}1Côté et al. (2001) adopted =1350 km/s., is the projected distance of each GC and is its position angle, measured East to North. We have used LevenbergMarquardt fitting method (Levenberg 1944; Marquardt 1963) to solve for and . 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.
4 Results and interpretation
4.1 Properties of globular clusters in the four groups
The cluster analysis divided the GCs of M87 into four
groups. We did not use the kinematical data ( and )
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 LevenbergMarquardt 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 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 Heenriched stars. So the interpretation of Fig.6, which
shows the run of H 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 /Fe. An
alternative explanation in terms of stellar populations with
younger ages is not possible because it would not be compatible
with the observed broadband colours (see below).
We next investigate the colours of the stellar
populations in the GCs of M87. The colourcolour 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 broadband 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 = 0.3 are shown as solid
lines, green and blue for 0 and 70% secondary populations
respectively. Changing 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 CFHMegcam
magnitudes to
SDSS using the paper by Betoule et al. (2013).
The colours of the metalpoor GCs (shown on Fig.7) are
well fit by the models. On the other hand, the metalintermediate
(G1) and metalrich (G2) GCs can only be explained by a
Heliumenriched 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
Heenriched 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 metalintermediate and
metalrich 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 superHerich 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 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 G1G3 and that the Na  O anticorrelation
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 (FraixBurnet et al. 2009), M31 (Colucci et al. 2012a) and the LMC (Colucci et al. 2012b), as expected from arguments of nucleosynthesis. High values of /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 metalpoor environments. In M87, we find that the opposite is true. We use Mgb/Fe = 2Mgb/(Fe5270 + Fe5335) as an indicator of /Fe and calibrate this indicator with the stellar population models of Thomas et al. (2011). As shown on Fig.9, the lightelement 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 metalpoor, the disagreement between /Fe and its indicator Mgb/Fe suggests that the lightelement 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 nonrotating 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.
4.2 Origin of the four groups
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 metalricher material
which ends up in circular orbits, whereas the metalpoor 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 metalpoor material, G4 in the inner part from metalrich
material. The lack of rotation in G3, the flattened spatial
distribution, strong rotation and high average /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, ). 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 : 374.93 km for 300, and 179.92 km) for 300.
5 Conclusion
We have performed a multivariate analysis of the
chemical and photometric properties of the GCs of M87, and found
four groups of clusters with distinct properties. The latter are
used to propose a formation scenario for the different groups, in
terms of dissipational collapse of the galaxy, followed by
dissipationless accretion of matter.
Regardless of the clustering results, the reanalysis of data from the literature with recent stellar population models, those of Thomas et al. (2011) and the Yonsei models (Chung et al. 2013), provides new insights into the chemical properties of the GCs of M87, compared to those of other galaxies. We have found that these properties are different from those of GCs in M31, in that the Lick indices NaD and TiO1 are higher at a given Mgb. We have also found a progressive excess blueing of the GCs with metallicity, which we interpret by the presence of an extreme blue HB, possibly due to He enriched secondary stellar populations. Finally, /Fe increases with metallicity, at variance with what occurs in GCs of spiral galaxies (our own and M31). Our interpretation is one among many others and remains speculative in view of the smallness of the sample and is not an ambitious as is needed. Our interpretation is only one possible way of understanding the results.
DA groups  G1  G2  G3  G4 

51  2  3  0  
2  52  1  1  
2  0  23  0  
0  0  0  10  
Total  55  54  27  11 
Parameter  G1  G2  G3  G4 

No of GCs  55  54  27  11 
2.02710.0735  1.94890.08  2.06070.08  1.2380.13  
Mg1  0.009070.004  0.012040.004  0.03230.005  0.10150.01 
Mg2  0.086180.007  0.092700.007  0.019810.009  0.27810.02 
Mgb  2.1230.115  2.5370.12  1.15110.09  4.7650.19 
Fe5270  2.0190.102  1.93810.073  0.7110.13  2.8760.17 
Fe5335  1.54820.09  1.64810.08  0.54480.09  2.4050.24 
NaD  1.51220.09  1.66460.001  1.030.11  5.1250.52 
TiO1  0.011690.002  0.01720.04  0.000810.002  0.036820.003 
Fe/H  0.98070.05  0.8320.07  1.62890.08  0.060.06 
Mgb/Fe  1.2318 0.0628  1.41840.0513  3.491.36  1.85430.0925 
20.2290.06  20.0770.006  20.4760.06  20.0020.13  
gr  0.602910.005  0.618330.01  0.550370.006  0.720.02 
gi  0.87270.01  0.90760.01  0.78370.01  1.08550.04 
147.738.36  338.038.92  234.217.9  146.419.6  
30.92 7.42  232.96 69.73  40.93 10.34  249.02 19.11  
92.58 14.27  265.38 19.41  167.20 13.59  217.14 50.65 
6 Acknowledgements
We thank C. Chung for sending us tables for the Yonsei Evolutionary Population Synthesis models. We are also thankful to A.K.Chattopadhyay for useful discussions, and to Saptarshi Mondal for helping in some calculations. One of the authors (Tanuka Chattopadhyay) thanks DST, India for providing her a major research project for the work.
References
 Albazzaz et al. (2004) Albazzaz, H., Wang, X.Z., 2004, Industrial & Engineering Chemistry Research, 43(21), 6731.
 Babu et al. (2009) Babu, G.J, Chattopadhyay,T., Chattopadhyay, A.K., Mondal, S., 2009, ApJ,700,1768.
 Bekki & Chiba (2005) Bekki, K., Chiba,M., 2005, ApJ, 625, L107.
 Betoule et al. (2013) Betoule,M., Marriner, J.; Regnault, N.; Cuillandre, J.C.; Astier, P.; Guy, J.; Balland, C.; et al. 2013, AA, 552, A124.
 Brodie & Strader (2006) Brodie,J.P., Strader,J., 2006, ARA&A 44, 193.
 Chattopadhyay et al. (2013) Chattopadhyay, T., Karmakar, P., 2013a, New Astron. 22, 22.
 Chattopadhyay et al. (2013) Chattopadhyay, A. K., Mondal, S., Chattopadhyay, T., 2013b, CSDA, 57, 17
 Chattopadhyay et al. (2009) Chattopadhyay, A.K., Chattopadhyay, T., Davoust, E., Mondal, S., 2009, ApJ, 705, 1533
 Chattopadhyay et al. (2010) Chattopadhyay, T., Sharina, M., Karmakar, P., 2010, ApJ, 724, 678
 Chattopadhyay et al. (2012) Chattopadhyay, T., Sharina, M., Davoust, E., De,T., Chattopadhyay, A.K., 2012, ApJ, 750, 91.
 Chattopadhyay et al. (2007) Chattopadhyay, T., Chattopadhyay, A. 2007, A&A, 472, 131
 Chung et al. (2013) Chung, C., Lee, S.Y., Yoon, S.K, Lee, Y.W, 2013, ApJ,769, L3
 Cohen & Ryzhov (1997) Cohen, J.G., Ryzhov, A., 1997, ApJ, 486, 230
 Cohen et al. (1998) Cohen, J.G., Blakeslee, J.P., Ryzhov.A., 1998, ApJ, 808, 539.
 Cohen et al. (2003) Cohen, J.G., Blakeslee, J.P., Côté, P., 2003, ApJ 592, 866
 Colucci et al. (2012) Colucci, J.E., Bernstein, R.A., Cameron, S.A., McWilliam, A., 2012, ApJ, 746, 29
 Colucci et al. (2012) Colucci, J.E., Bernstein, R.A., Cohen, J, 2012, Proc. XII International Symp. on Nuclei in the Cosmos, 2012, 99
 Comon (1994) Comon, P., 1994, Signal Processing, 36, 287.
 Cote et al. (2001) Côté, P., McLaughlin, D.E., et al, 2001, ApJ, 559, 828.
 Eggen et al. (1962) Eggen, O.J., LyndenBell, D., Sandage, A. R., 1962, ApJ 136, 748
 Emselhem et al. (2007) Emselhem, E., Cappellari, M., Krajnovic, D., et al, 2007, MNRAS, 379, 401
 FraixBurnet et al. (2010) FraixBurnet, D., Chattopadhyay, T., Chattopadhyay, A.K., Davoust, E., 2010, MNRAS, 407, 2207
 FraixBurnet et al. (2010) FraixBurnet, D., Davoust, E. and Charbonnel, C. 2009, MNRAS, 398, 1706.
 Hanes et al. (2001) Hanes, D.A., Cöté, P., Bridges, T.J., McLaughlin, D.E., Geisler, D., Harris, G.L.H., Hesser, J.E., Lee, M.G., 2001, ApJ 559, 812
 Huchra & Brodie (1987) Huchra, J.P., Brodie, J.P., 1987, AJ, 93, 779.
 Jesseit et al. (2009) Jesseit, R., Cappellari, M., Naab, T., et al, 2009, MNRAS, 397, 1202.
 Johnson & Wichern (1998) Johnson, R. A., & Wichern, D. W. 1998, Applied Multivariate Statistical Analysis (4th ed.; Upper Saddle River: Prentice Hall)
 Kaviraj et al. (2007) Kaviraj, S.; Sohn, S. T., O’Connell, R. W., Yoon, S.J., Lee, Y. W., Yi, S. K., 2007, MNRAS, 377, 987
 Kotulla et al. (2009) Kotulla, R., Fritze, U., Weilbacher, P., & Anders, P., 2009, MNRAS 396, 462
 Krajnovic et al. (2008) Krajnovic, D., Bacon, R., Cappellari, M., Davies, R. L., de Zeeuw, P. T., Emsellem, E., FalcónBarroso, J., et al. 2008, MNRAS, 390. 93.
 KisslerPatig & Gebhardt19 (1998) KisslerPatig, M. & Gebhardt, K. 1998, AJ, 116, 2237.
 Lada & Lada (2003) Lada, C.J., Lada, E.A., 2003, ARA&A, 41, 57.
 Lambert & Gontier (2009) Lambert, S.B., Gontier, A.M., 2009, A&A 493, 317.
 Levenberg (1944) Levenberg, K., 1944, Q Appl. Math., 2, 164.
 Macqueen (1967) Macqueen, J., 1967, in Proc 5th Berkeley Symp. Math Prob. Vol.1, ed. L.M. Lecom and J. Neyman (Los Angeles Univ of California press), 281.
 Maratson et al. (2003) Maratson, C., Greggiol, Renzini, A., et al, 2003, AA, 400, 823.
 Marquardt (1963) Marquardt, D.W., 1963, SIAMJ, Appl. Math., 11, 431.
 Martell & Grebel (2010) Martell, S.L., and Grebel, E.K., 2010, AA, 519, 14.
 Montes (2014b) Montes, M., Trujillo, I., et al, 2014b, MNRAS, 439, 990.
 Montes et al. (2014a) Montes, M., Trujillo, I., Prieto, M.A., AcostaPulido, J.A., 2014a, MNRAS, 439, 990.
 Muratov et al. (2010) Muratov, A.L. Gnedin, O.Y., 2010, ApJ, 718, 1266.
 Mould et al. (1987) Mould, J.R., Oke, J.B., Nemec, J.M., 1987, AJ, 93, 53.
 Mould et al. (1990) Mould, J.R., Oke, J.B., De Zeeuw, P.T., Nemec, J.M., 1990, AJ, 99, 1823.
 Park et al. (2012) Park, H.S., Lee, M.G., Hwang, H.S., Kim, S.C., Arimoto, N., Yamada, Y., Tamura, N., Onodera, M., 2012, ApJ 759, 116
 Proctor et al. (2009) Proctor, R.N., Forbes, D.A., Romanwosky, A.J., et al, 2009, MNRAS, 398, 91.
 Richtler et al. (2004) Richtler, T., Dirsch, B., Gebhradt, K., Geisler, D., et al, 2004, AJ, 127, 2094.
 Schubert et al. (2010) Schuberth Y., Richtler T., Hilker M., Dirsch B., Bassino L.P., Romanowsky A.J., Infante L., A&A 2010, 513, A52
 Schiavon et al. (2012) Schiavon, R. P., Caldwell, N., Morrison, H., Harding, P., Courteau, S., MacArthur, L. A., Graves, G. J. 2012, AJ, 143, 14.
 Shapiro et al. (1965) Shapiro, S.S., Wilk, M.B., 1965, Biometrika, 52 , 591.
 Smith et al. (2000) Smith, R. J., Lucey, J. R., Hudson, M. J., Schlegel, D. J., Davies, R. L. 2000, MNRAS, 313, 469
 Sohn et al. (2006) Sohn, S. T., O’Connell, R. W., Kundu, A., Landsman, W. B., Burstein, D., Bohlin, R. C., Frogel, J. A., Rose, J. A., 2006, AJ 131, 866
 Strom et al. (1981) Strom, S.E., Strom, K.M., Wells, D.C., Forte, J.C., Smith, M.G., Harris, W.E., 1981, ApJ, 245(5457), 416.
 Strader et al. (2011) Strader, J., Romanowsky, A. J., Brodie, J. P., Spitler, L. R., Beasley, M. A.; Arnold, J. A., Tamura, N.et al. 2011, ApJS, 197, 33.
 Sugar & James (2003) Sugar, A.S., James, G.M., 2003, J.Am. Stat. Assoc., 98, 750.
 Thomas et al. (2011) Thomas D., Maraston C., Johansson, J., 2011, MNRAS 412, 2183
 Woodley et al. (2007(@) Woodley, K.A., Harris, W.E., Beasley, M.A., et al, 2007, AJ, 134, 494.