Chemodynamical Clustering Applied to APOGEE Data: Rediscovering Globular Clusters

Globular clusters are characterized by stars with significant anticorrelations among [Al/Fe]–[Mg/Fe] and [Na/Fe]–[O/Fe] but overall homogeneous Fe (Carretta et al. 2009a, 2009b). Before describing the algorithms that perform the clustering analysis on the data set, we will investigate the presence of these anticorrelations in the DR12 Cannon data for the globular clusters M107, M13, M15, M2, M3, M5, M53, M92, and NGC 5466.

Figure 1 plots [Na/Fe]–[O/Fe] (top panels) and [Al/Fe]–[Mg/Fe](bottom panels) for stars in the DR12 Cannon that are members of the observed globular clusters. M107, M5, M13, M3, and M2 show evidence of anticorrelation in the [Al/Fe]–[Mg/Fe] plots. However, [Na/Fe] does not anticorrelate with [O/Fe] in the DR12 Cannon, except for a weak anticorrelation displayed in M92. One possible reason for the lack of anticorrelation between Na and O is that Na has relatively few and weak features, as discussed in Holtzman et al. (2015). Fernández-Trincado et al. (2017) reported peculiar abundance patterns in 11 red giant stars that show Al and N enhancements with C and Mg depletions. These patterns are similar to those exhibited in second-generation stars of globular clusters. Similarly, Martell & Grebel (2010) and Martell et al. (2011) found globular-cluster-like C–N abundance variations in a subset of halo field stars, suggesting that these stars had a globular cluster origin.

Zoom In Zoom Out Reset image size

In Figure 1, we also found an anticorrelation between [Al/Fe] and [Mg/Fe] in some background stars with [Fe/H] < −1. These background stars correspond to members of the M13 and possibly M2, M3, and M5 globular clusters in the chemical space. Ideally, a successful chemical tagging algorithm should be able to recover this overdense structure. However, as we noted before, our results remained largely unchanged when all elements were included in the algorithms. Even though these background stars appear as an overdense region when projected onto [Al/Fe] and [Mg/Fe], conditioned on their [Fe/H] values, our ability to recover them with chemical tagging is limited. These stars do not share similar radial velocities and can potentially be separated with chemodynamical tagging.

M5, M15, and M107 are the only globular clusters that are observed by APOGEE and also present in the sample studied by Carretta et al. (2009a) at optical wavelengths. When comparing our [Al/Fe]–[Mg/Fe] plots with those of Carretta et al. (2009a), we noticed that the DR12 Cannon [Mg/Fe] abundances are systematically lower by ∼0.2 dex on average than those in Carretta et al. (2009a), while the ranges in [Al/Fe] are consistent with each other. Moreover, while the data collected by Carretta et al. (2009a) for red giant stars in these clusters do not show any kind of correlation, we see a clear [Al/Fe]–[Mg/Fe] anticorrelation in M5 and possibly in M107, but nothing in M15 in our plots.

We also compared the DR12 Cannon distributions with those obtained by Mészáros et al. (2015). Both data sets cover the same abundance ranges for each cluster except for NGC 5466, whose DR12 Cannon [Mg/Fe] abundances are shifted toward lower values by ∼0.15 dex. Both Carretta et al. (2009a) and Mészáros et al. (2015) confirm an [Al/Fe]–[Mg/Fe] anticorrelation in globular clusters. The Cannon data indicate a possible [Al/Fe]–[Mg/Fe] anticorrelation in M107, which is not seen in Mészáros et al. (2015), but do not exhibit any [Al/Fe]–[Mg/Fe] anticorrelation in NGC 5466, M15, M53, and M92, otherwise detected in Mészáros et al. (2015).

As a first step, we used the K-means algorithm in order to rebuild globular clusters that are expected to appear as overdense or clustered regions in the N-dimensional chemical space. The K-means algorithm is a clustering algorithm in data mining that has been applied for chemical tagging in previous work (Hogg et al. 2016). The algorithm aims to partition data points into K groups, each of which contains a core position (not necessarily a data point) and the data points nearest to it.

Here we explain how K-means operates in the N-dimensional abundance spaces following Hogg et al. (2016):

• 1.  
  Choose a value of K that defines the number of groups in the N-dimensional space.
• 2.  
  An initial position is chosen to be the center for each group.
• 3.  
  Assign each star to the group with the closest center.
• 4.  
  Recenter each group to the centroid of its stars.
• 5.  
  Iterate the procedure in steps 3 and 4 until convergence has been reached.

As discussed before, we chose to work in the iron-normalized abundance space. Compared to Hogg et al. (2016), we discarded six abundances according to anticorrelations found in [Al/Fe]–[Mg/Fe] and [Na/Fe]–[O/Fe], discussed in Carretta et al. (2009a), as well as [C/Fe] and [N/Fe]. K-means is usually run several times, with different K to determine the optimal number K. However, we only show the results with K = 512 here. Other choices of K do not significantly improve the results.

While widely popular, K-means faces some limitations. In its initial and basic implementation, it uses Euclidean distances in the N-dimensional space and assumes a spherical shape and a fixed size of potential clusters in the data. In particular, the K-means algorithm assigns every star to one group. This limits its applicability, since background noise stars are forced by the algorithm to always be assigned to a group. This limitation also means that the overdensities recovered by K-means in the chemical space will have many false positives and are not guaranteed to be true clusters. We find that for low values of K the algorithm tends to combine separate known clusters into a mono-abundance group. At high K values, K-means tends to divide a known cluster into separate groups.

Figure 2 displays the two-dimensional projections of our sample in the chemical space and the results from the K-means algorithm with K = 512. Colored symbols represent the APOGEE stars that are known members of globular clusters (with different colors representing the different clusters targeted in the data set). The gray symbols are the background stars (not members of any globular cluster). The black encircled symbols are the members recovered by K-means. This figure demonstrates that stars appearing to be similar in chemical properties according to the K-means algorithm are rarely members of any known globular clusters. Members of the M13 globular cluster are recovered by the algorithm with a recovery rate of 26.7% and a false-positive rate of 66.7% in a single K-means group. None of the members in the other globular clusters are identified in our sample with K-means.

For a general data set, if we know the number of disjoint groups and that data points within each group are normally distributed in our data, the underlying structures of the data set can usually be determined by using K-means. However, in many cases the groups that characterize the data set do not have a spherical shape, or their number is unknown. In these cases, a more flexible algorithm is required to rebuild globular clusters (e.g., DBSCAN), to account for more complex shapes and sizes of clusters. The basic premise behind the algorithm is that each group is expected to be defined by a probability distribution and the probability of picking a star near the central core of the distribution should be much higher than at the edge. If two groups overlap slightly, the overlap will be within the outer edges of the distributions, where the probability is negligibly low. Hence, DBSCAN ignores the data points that are sparse and focuses on connecting the dense regions, which should be the central cores of potential groups.

Two parameters have to be chosen for this algorithm: the maximum radius of the neighborhood from a core position, , which determines how large an area around each data point is considered; and the minimum number of data points K inside the neighborhood for a point to be considered a core point, MinPts. A data point is a core point if there are at least K other data points that are at most a distance away from it. The DBSCAN algorithm works by connecting all pairs of data points such that the two data points (stars) are at most away and at least one of them is a core point.

This algorithm has some limitations as well. The results rely on the choice of parameters and K. The results can be improved by running the algorithm multiple times on data with different parameters. The more important caveat is that it does not perform effectively on data sets where the overlap of the two groups extends more than just along their periphery, or where distinct groups have different internal densities. In the latter case, DBSCAN might interpret the less dense group as a collection of outliers. However, the definition of the periphery changes with different norms. Therefore, it is possible to overcome this limitation by adopting a different norm. Consistent with Mitschang et al. (2013), we use the Manhattan distance, or L₁ norm, as the metric for implementing chemical tagging with DBSCAN.

Figure 3 illustrates the two-dimensional projections of the chemical space for our sample and the recovered members according to the DBSCAN algorithm with  = 0.4 and K = 6. We found that the recovery rate quickly diminishes when we deviate from these values. As compared to chemical tagging performed with K-means, DBSCAN recovers M13 at a similar rate (23.9%) and 18.4% of M5 members. However, the M13 members are recovered in four DBSCAN groups, compared to just one group for K-means.

Previous sections showed that the process of identifying groups in chemical space can retrieve at most 30% of true members in known globular clusters. This result holds independently of the algorithm adopted to identify clusters in the chemical space. This is likely due to insufficient precisions of chemical abundances. We found that radial velocities of members from distinct globular clusters are almost indistinguishable in several cases. Thus, radial velocity alone is likely not sufficient to recover individual globular clusters either.

As the next step, we add the kinematic information, i.e., radial velocity, to the chemical space and propose a new algorithm, SNN. We used AstroPy (Astropy Collaboration et al. 2013) and various SciPy packages (Jones et al. 2001), including Numpy (Walt et al. 2011) and Scikit-learn (Pedregosa et al. 2011), to build our algorithm. We applied this algorithm to our sample and identified star clusters in the chemistry–velocity space.

All clustering techniques rely on assigning stars based on some distance metric. One problem with picking a distance metric is normalization when two spaces with different ranges are considered (e.g., the chemical space and radial velocity space). One of the main advantages of our algorithm is that it identifies clusters in the chemistry–velocity space by taking the intersection of the nearest neighbors in the chemical space and radial velocity space. This feature avoids normalization between chemical abundances and radial velocity and guarantees that the final nearest neighbors are neighbors in both the chemical and velocity space.

The algorithm operates as follows. As suggested in Mitschang et al. (2013), for each set of elements $\{{X}_{i}| i=1,2,\,\ldots ,\,N\}$, the distance in the chemical space between two stars, A and B, is computed as the Manhattan distance, defined below:

$$\begin{eqnarray}&&{d}_{\mathrm{chem}}=\displaystyle \sum _{i=1}^{N}| {X}_{{iA}}-{X}_{{iB}}| ,\end{eqnarray} \tag{ 1 }$$

where X_i = [Fe/H], [K/Fe], [S/Fe], [Si/Fe], [Ca/Fe], [Mn/Fe]...

The distance between the radial velocities of two stars is a simple Euclidean distance:

$$\begin{eqnarray}&&{d}_{\mathrm{rv}}=\sqrt{{({v}_{A}-{v}_{B})}^{2}}.\end{eqnarray} \tag{ 2 }$$

For every star, a list of nearest neighbors in the chemical space and a list of nearest neighbors in the radial velocity space are respectively found. The neighbors for a particular star in the chemodynamical space are the intersection of neighbors in the chemical space and velocity space.

Then, the Jaccard distance between each pair of stars is computed by

$$\begin{eqnarray}&&{d}_{J}({G}_{A},{G}_{B})=1-\displaystyle \frac{| {G}_{A}\cap {G}_{B}| }{\left|{G}_{A},\cup ,{G}_{B}\right|},\end{eqnarray} \tag{ 3 }$$

where G_A and G_B are the lists of neighbors for star A and B, respectively, ∩ is the intersection, ∪ is the union of two sets of neighbors, and the absolute values represent the cardinality of a set. The chemistry–velocity clustering analysis is performed next using the Jaccard distance values and the DBSCAN algorithm.

Our approach has two major advantages. First, as mentioned earlier, we are able to utilize radial velocity and chemical abundances together without compromising either. Each neighbor of a star is a neighbor in both chemical and velocity space. Second, unlike K-means and DBSCAN, the SNN algorithm does not make assumptions about the sizes, shapes, and densities of clusters (Ertoz et al. 2003).

A schematic view of our methodology is illustrated in Figure 4 and summarized as follows:

• 1.  
  Calculate for each pair of stars the distance d_chem in the N-dimensional chemical space (N = 9).
• 2.  
  Calculate for each pair of stars the distance d_rv in radial velocity.
• 3.  
  Find the N_chem nearest-neighbor stars in the chemical space.
• 4.  
  Find the N_rv nearest-neighbor stars in the radial velocity space.
• 5.  
  Take the intersection of neighbors in the chemical space and velocity space.
• 6.  
  Remove stars with fewer than N_cut neighbors.
• 7.  
  Calculate the Jaccard distance between every pair of stars.
• 8.  
  Perform DBSCAN clustering with the precomputed Jaccard distances.

Zoom In Zoom Out Reset image size

Based on exploratory analysis, we chose N_chem = 700 and N_rv = 500. We picked both numbers to ensure that most known members of globular clusters have some neighbors in the chemistry–velocity space. Given the choices of N_chem and N_rv, we found that 96% of globular cluster members have at least one neighbor in the chemistry–velocity space. We will not have this information in a true blind-search scenario. We instead have to rely on other methods to choose parameters. The numbers N_chem and N_rv should also not be set too high in order not to cause members from separate globular clusters to share neighbors. We choose a higher number for N_chem because globular clusters in the DR12 Cannon are more homogeneous in radial velocity than in chemical abundances. The parameters used for DBSCAN are  = 0.35 and MinPts = 7. These parameters indicate that a core point should have at least seven neighbors that share at least 65% common neighbors out of all neighbors. If a star has fewer than N_cut neighbors, it will be removed from our sample before we perform DBSCAN. We used N_cut to eliminate stars without any neighbors in our implementation, or equivalently N_cut = 1.

We show next that our methodology based on the SNN algorithm is promising in discovering new globular clusters or detecting debris stars related to the globular clusters orbiting in our Galaxy.

We applied the SNN algorithm to the DR12 Cannon sample, searching for clusters in the chemistry–velocity space. The groups identified with the SNN are then matched in equatorial coordinates with the catalog of Galactic globular clusters of Harris (1996). The following globular clusters have been identified: M13, M3, M92, and NGC 5466, for which we retrieved the corresponding tidal radii as listed in Harris (1996). However, the membership information only comes from APOGEE, consistent with our previous discussion. Similar to Figures 2 and 3, we plot the recovered stars from SNN in Figure 5. Although we highlight members recovered in groups with a false-positive rate lower than 85%, the exact false-positive rates are shown in Table 1.

Zoom In Zoom Out Reset image size

Our implementation of the SNN algorithm took ∼20 minutes. Although it is slower than DBSCAN alone (∼5 minutes), it is significantly faster than the standard implementation of K-means in scikit-learn (∼75 minutes).

Once the groupings were returned with their stars identified as candidate members, we displayed for each group the two-dimensional projections in the chemical space. An example is illustrated by group 32 identified by the SNN algorithm as a cluster of stars in the chemistry–velocity space. Although it does not overlap perfectly with a globular cluster, SNN group 32 is dominated by members of M13, which are indicated by open circles in Figure 5.

Figure 6 shows that stars of group 32 display an anticorrelation between [Al/Fe] and [Mg/Fe], which is a clear characteristic of globular clusters, even though our algorithm did not use [Al/Fe] and [Mg/Al]. A subsequent inspection shows that this clump of stars is indeed the globular cluster known as M13, found without a priori knowledge of its location in the sky. Similar to Figure 5, the stars recovered by the algorithm and confirmed to be members are circled with open symbols in Figure 6. Note that our algorithm identifies more stars than the known confirmed members.

Zoom In Zoom Out Reset image size

A study of these nonconfirmed members shows evidence of anticorrelation between [Al/Fe] and [Mg/Fe]. It should be noticed, however, that while the majority of the nonconfirmed members are within one tidal radius of the globular cluster (orange points), the remaining stars do follow the anticorrelation but are located outside the tidal radius. Two stars are between one and two tidal radii (blue triangles), and 15 stars are beyond three tidal radii away (green symbols). It is possible that these stars are candidate stellar debris of the globular cluster. Interestingly, even though our algorithm did not include [Al/Fe] and [Mg/Fe], we recovered the overdense region in [Al/Fe] and [Mg/Fe] mentioned in Figure 1. Conversely, when we do include [Al/Fe] and [Mg/Fe] in our chemical tagging algorithms, our recovery rates never exceed 30% for M13, much lower than the 95.8% recovery rate from chemodynamical tagging.

As a next step, we cross-matched the sample of stars obtained with our chemodynamical tagging algorithm with the Two Micron All Sky Survey (2MASS) catalog (Skrutskie et al. 2006). Figure 7 displays the resulting photometry on a color–magnitude diagram (CMD) together with the PARSEC isochrone (Bressan et al. 2012) of M13 calculated with the parameters of O'Malley et al. (2017) (E(B − V) = 0.02, [Fe/H] = −1.58, (m − M)_V =14.54 mag, and age = 12.3 Gyr). The symbols and colors are the same as in Figure 6. In the bottom right corner of the plot we show the mean photometric errors of the sample. The distribution of the stars within one tidal radius on the CMD is in good agreement with the isochrone representing the cluster, and the stars with distances up to three tidal radii also fall within the same distribution. In addition, 5 out of 15 stars at distances larger than three tidal radii are consistent with the infrared photometry of M13. While most outliers lie on the red side of the CMD, two stars are bluer than the mean (J–K_s) color of the cluster, and one of them is situated within one tidal radius, 2M 16414108+3627550. The unexpected color of this star is explained in the 2MASS archive: the J-band flux for this star is flagged as "detected at the location where the images were not deblended consistently in all three bands (JHK_s)" and does not have any error listed in the 2MASS catalog. Unfortunately, we could not perform a photometric analysis of the data at visible wavelength because the Hubble Space Telescope observations of M13, although available, constrain themselves at the very center of the cluster and thus exclude part of our sample.

Zoom In Zoom Out Reset image size

We recovered 34.2% of M3 members, with a false-positive rate of 34.9%, and 75% of NGC 5466 members, with a false-positive rate of 14.3%, each in their respective groups.

We did not recover M5, even though it has 103 members in our data set. There are too many background stars with similar radial velocities to members of M5. There are 13,878 stars within the range of radial velocities of M5 members. Our choice of N_rv = 500 is likely too low for M5 members to be the nearest neighbors of each other in the velocity space and, in turn, chemistry–velocity space.

We recovered known members from the M92 cluster, but they were split into two separate groups. The two groups contain 26 and 10 members, respectively, with false-positive rates of 44.7% and 44.4%. In Table 1, we show the recovery and false-positive rates for the two groups combined. The split of M92 could be caused by the large number of background stars with similar radial velocities. There are 2063 stars within the range of radial velocities of the M92 members. Again, our choice of N_rv might be too low for the M92 members to be connect in velocity space.

When N_rv is increased, we notice that different globular clusters start to merge with each other in our results. Since radial velocity is only one-dimensional, the data points are not able to spread out, as in the N-dimensional chemical space. We expect this to change when radial velocity is replaced with a higher-dimensional velocity space.