Tracing Birth Properties of Stars with Abundance Clustering

We focus our analysis on the present-day disk. We first select stars that overlap in space with the disk by imposing limits of ∣z∣ ≤ 0.5 kpc and 4 kpc ≤ R_GAL ≤ 12 kpc, though our results are consistent for other spatial cuts. We then determine stars that are current disk members using three-dimensional velocity space. We model a two-component Gaussian mixture model in (v_θ , v_r , v_z ), similar to the approach taken by Buck et al. (2019b) and Obreja et al. (2018a) to model simulations in kinematic space using Galactic Structure Finder (GSF; Obreja et al. 2018b), and define disk stars to have a Mahalanobis distance less than 2 from the center of the corresponding Gaussian distribution. Figure 1 shows these selection cuts in the equivalent Toomre diagram and x − y plane for both simulations. Finally, we do an additional cut of R_birth ≤ 15 kpc to ensure we are not looking at infalling debris.

Zoom In Zoom Out Reset image size

Additionally, since the abundances cover different ranges, we make quality cuts on our abundance data in the scaled ([Fe/H], [X/Fe]) space, where the transformed abundances have mean 0 and a standard deviation of 1. Since the goal of this paper is to focus on global properties between abundances, age, and R_birth, we remove outliers by only selecting stars that have scaled abundances between –4 and 4. Our final sample sizes are 229,045 and 44,359 particles for Chem2d and Chem15d, respectively.

Figure 2 shows the simulation data in the [α/Fe]–[Fe/H] plane after the selection cuts discussed above. Due to the formation history, both simulations have obvious trends in R_birth and age (middle and right columns). For a given value of [α/Fe], R_birth increases as [Fe/H] decreases. Similarly, for a given value of [Fe/H], age increases as [α/Fe] increases. As discussed in Buck (2020), the horizontal age gradient and the diagonal radius separation in [α/Fe]–[Fe/H] are simply a reflection of star formation in the disk happening at different radii, where metallicity decreases with increasing radii.

Zoom In Zoom Out Reset image size

The left column of Figure 2 shows the density structure of the two simulations. For Chem2d (top-left), there are density ridges which follow along different age tracks, whereas there is no noticeable substructure in Chem15d (bottom-left), presumably due to its lower mass resolution. The left panel of Figure 2 also shows that the footprint of Chem15d is different than Chem2d. Most noticeably, the spread in this plane primarily captures high-α stars for Chem15d. The structural differences between the [α/Fe]–[Fe/H] planes of Chem2d and Chem15d are due to slight differences in their formation history (discussed above) and the different set of chemical yields for chemical enrichment (see Buck et al. 2021 for discussion on the impact of yield tables and tracks in this abundance plane).

Following the same methodology Ratcliffe et al. (2020) used with observational red clump DR14 apogee data, we use agglomerative hierarchical clustering using Ward's minimum variance criterion (Ward 1963) as one of the ways to combine the most chemically similar stars. Specifically, we use the Ward2 algorithm described in Kaufman & Rousseeuw (2009) and Murtagh & Legendre (2014). We conceptually describe the algorithm here, and refer the reader to Ratcliffe et al. (2020) for a more in-depth explanation.

The algorithm begins with each star as its own cluster, and at each step we combine the pair of clusters that leads to a minimum increase in total within-cluster variance until only one large cluster containing all the stars remains. The output is a tree showing the combination of groups at each step, called a dendrogram. Thus, the user decides the number of clusters to separate the data into after seeing the linking structure of the data.

EnLink is a nonparametric hierarchical clustering algorithm built on a locally adaptive distance metric and hence able to identify complex structures in the data. The full data set is first divided via a binary-partitioning algorithm which uses an entropy criterion to preferentially bisect dimensions that contain maximum information ("EnBid"; Sharma & Steinmetz 2006). This approach allows a nonparametric definition of "local" regions in the data set from which the adaptive metric, with flexible scales and orientations, can then be derived (see Sharma & Johnston 2009 for full details). It is this metric that defines the distance between particles subsequently.

EnLink partners the adaptive distance metric machinery with IsoDen (Pfitzner et al. 1997), which is a density-based clustering algorithm. Conceptually, clusters can be thought of as regions around high-density peaks that are separated from one another by lower-density regions. Thus, as we lower the isodensity contours when examining a high-density region, we stay within the cluster until we encounter a lower-density region that connects to another high-density region. Then, as the isodensity contour continues to lower, a new group encompassing both clusters is formed. Continuing in this fashion forms a hierarchy of density-based parent–child clusters.

EnLink has two user-specified parameters: the number of nearest neighbors used in calculating density (k_den) and the threshold significance level when comparing the high and low density levels of parent–child clusters (S_th). Since the goal of our analysis is not to find the best clustering but rather to investigate the stability of our results, we choose to vary k_den from 30 to 1000, and S_th is such that the expected number of groups due to Poisson noise is 0.5, 1, and 2. We did not observe major differences in our results.

As mentioned in Section 3.1, hierarchical clustering produces a dendrogram showing how stars in the abundance space combine, starting from each star being its own group to one cluster containing every star. After the linking structure is determined, the user then specifies how many groups to separate the sample into. Walking down the tree, and thus increasing the number of clusters, corresponds to one group being separated to form two new clusters at each step. For a given number of k groups defined in the two- or 15-dimensional abundance space, we determine the contour level that contains 50% of the stars within each abundance group after projecting into the age–R_birth plane. We then calculate the percent that each 50% contour level overlaps with the other k − 1 group's 50% contour levels by laying down a fine grid and comparing the number of points in just the ith group to the number of points that fall in more than just the ith group. For the ith group, the overlap percentage is defined as the percent of area that is common between the 50% contour region of the ith group and the 50% contour regions of the other k − 1 groups.

The top-left panel of Figure 3 shows the median of these overlap percentages and standard deviation of the k overlap percentages as a function of the number of groups found in Chem2d using hierarchical clustering. We see that groups defined solely using two abundances show consistently low overlap in birth time and space for up to 10 groups at the 50% contour level (and up to seven groups at the 75% contour level, which is not shown).

Zoom In Zoom Out Reset image size

The middle and right panels of the top row in Figure 3 show the seven groups in the [α/Fe]–[Fe/H] and R_birth–age planes at the 50% contour level. We can see that the groups found using two abundances separate diagonally, both as a function of age and R_birth. We believe that this primarily is a consequence of formation history, as the low-α stars have gradients in abundances, age, and R_birth.

4.1.1. Using EnLink to Leverage Density Structure in Chem2d

In the two-dimensional abundance space of Chem2d, we can see streaks of higher-density regions along the different age bins (see top-left of Figure 2). Hierarchical clustering does not leverage the density of the simulation in abundance space and thus is unable to capture the streak formations (see top-middle of Figure 3). Therefore, to include this structure as information in assignment of cluster groups, we employ a density-based clustering method with an adaptive distance metric to use the ridge-like structure in association of groups.

As discussed in Section 3.2, EnLink has the number of nearest neighbors used to determine the density at a point and the maximum number of spurious clusters created by noise as input parameters. We find that the group separation in the age–R_birth plane is fairly stable when we focus on allowing either 0.5 or 1 spurious clusters and 30 to 1000 nearest neighbors, with the majority of EnLink groups having a median overlap percent of below 25% at the 75% contour level, and near 0% at the 50% contour level.

Figure 4 shows the 50% contours in the abundance and age–R_birth planes for seven groups and parameter settings 260 nearest neighbors and a maximum of one spurious cluster. We can clearly see that EnLink captures the streaks of overdensity in the [α/Fe]–[Fe/H] plane better than hierarchical clustering, and in doing so we see better separation in the age–R_birth plane, with a median overlap percentage of only 1%. Since EnLink traces the higher-density ridges, and the ridges trace different age bins, the group separation in the age–R_birth plane no longer follows the diagonal trends given by hierarchical clustering. The right panel of Figure 4 shows that groups which follow along the same age track (e.g., the middle-aged blue, orange, and green groups), have similar ages and are separated as a function of birth radius. This causes the blue and green groups to flip from the groups found using hierarchical clustering, but this is an inconsequential difference between the two clustering methods. For an easier way to directly compare the results of hierarchical clustering and EnLink in Chem2d, please see Figure 18 in the Appendix.

Zoom In Zoom Out Reset image size

Each EnLink group primarily lives in a unique birth time and place, though there is some overlap with the middle-aged groups, possibly due to not fine tuning the algorithm. Overall, we conclude that for the high-resolution Chem2d simulation, using a density-based method with an adaptive metric is desirable for the best results.

So far we have demonstrated that just two abundances ([α/Fe] and [Fe/H]) can trace separate ages and birth radii. We now investigate how additional abundances, including additional nucleosynthetic families, strengthen this result. The list of abundances and their families is given in Section 2. Due to the difficult problem of estimating density in a high-dimensional space, in addition to the issue of tuning the algorithm, we do not use EnLink to cluster in 15 dimensions, and instead focus on hierarchical clustering.

As shown in the middle-left panel of Figure 3, groups defined in the chemical space of 15 abundances show more separation in the age–R_birth plane than the groups defined in the two-dimensional abundance space of Chem2d. These exhibit separation for 13 groups at the 50% contour level, and eight groups at the 75% contour level (not shown). As shown in the middle and right columns of the middle row of Figure 3, the groups comprised of older stars (which trace the high-α stars) primarily show separation as a function of age, whereas the younger low-α stars show separation in both R_birth and age. This shows that high-α stars were all born near the Galactic center, whereas the low-α stars were born at different radii and times throughout the galaxy.

Comparing these results to Chem2d (top-left of Figure 3), we see that groups defined in the 15-dimensional abundance space of Chem15d using hierarchical clustering overlap less in age and R_birth. In particular, for up to nine groups, the 15-dimensional groups are predominantly distinct in birth time and space at the 50% contour level, whereas in two dimensions the groups have some overlap even for as few as four groups.

In the previous section (Section 4.2) we showed that groups defined in the 15-dimensional abundance space of Chem15d showed more separation in age and R_birth than groups defined in the two-dimensional abundance space of Chem2d. Here we show the results of clustering in just the [α/Fe]–[Fe/H] plane of Chem15d are consistent to those of Chem2d.

The bottom-left of Figure 3 shows that the amount that abundance groups overlap in the age–R_birth plane is similar between Chem2d (red) and the [α/Fe]–[Fe/H] plane of Chem15d (black, labeled "projected"). As the simulations are split into more groups, both consistently show less distinction in (age, R_birth). On the other hand, the groups defined in the full 15-dimensional abundance space of Chem15d retain separate ages and birth radii for up to nine groups. This shows that more separate birth information is retained in abundance groups when more abundances are included, and therefore we claim abundance groups being more separate in age and R_birth is due to additional abundance information and not an artifact of different simulation history or resolution.

Figure 5 compares eight groups defined in the [α/Fe]–[Fe/H] plane of Chem15d to eight groups defined in the full 15-dimensional abundance space of Chem15d. The groups are arranged in order of age, with group 1 being the oldest and group 8 being the youngest. The oldest groups share the most stars between the two simulations, whereas the middle-aged and youngest groups are more muddled.

Zoom In Zoom Out Reset image size

So far we focused on how chemically similar stars differ in birth time and space using clustering methods. Now we give motivation for why use of a clustering method is needed for grouping stars.

Figure 6 shows the median and standard deviation of the percent group overlap for both Chem2d and Chem15d when separating the stars using different grouping methods. We compare hierarchical clustering with laying down a Cartesian grid in the [α/Fe]–[Fe/H] plane. The number of [α/Fe]–[Fe/H] bins are chosen to produce four, seven/eight (for Chem2d/Chem15d respectively), and 12 groups. We also show the results of EnLink for Chem2d. The left panel of Figure 6 shows that when only two dimensions of abundance information are known, we do not gain any more knowledge of birth properties from using a simple clustering method than if we were to create bins by laying down a Cartesian grid across the [α/Fe]–[Fe/H] plane. The groups retain more separate birth properties when leveraging density with an adaptive distance metric, however the errors are higher than that of binning and hierarchical clustering.

Zoom In Zoom Out Reset image size

The right panel of Figure 6 reveals that when higher-dimensional abundance information is available, there is a noticeable difference between only looking at stars separated using a grid in [α/Fe]–[Fe/H] versus using a clustering method in the full abundance space. Not only does clustering in 15 dimensions produce less overlap in the age–R_birth plane, but based on the smaller standard deviation, groups defined using hierarchical clustering produce consistently small overlap between groups whereas groups defined by binning in the [α/Fe]–[Fe/H] plane produce a large range of amount of groups that overlap in birth time and place.

Figure 7 compares the mean and standard deviation of (age, R_birth) of each group determined using two different methods: hierarchical clustering (left) and binning in the [α/Fe]–[Fe/H] plane (right) to create a Cartesian grid. This figure demonstrates that groups defined using hierarchical clustering preserve physical interpretation and agree with expected trends, namely we see a clear metallicity gradient as a function of age for a given R_birth when stars are separated using hierarchical clustering. The trend is not as obvious to see when stars are separated in the [α/Fe]–[Fe/H] plane with a grid. Additionally, we see that groups defined using a clustering method tend to have less overlap, and the overlap they do have is consistent among nearly all the groups. However, for the groups defined via a grid, the overlap between groups is irregular, with some groups being mainly separate and others completely overlapping multiple groups.

Zoom In Zoom Out Reset image size

Additionally, particularly for Chem15d , the age dispersion for each group when defined using hierarchical clustering is much smaller than when groups are separated using a grid. This again shows that the groups found using hierarchical clustering represent different physical groups in time (i.e., the groups have differing properties in birth age and location). Establishing a connection between those groups and star formation episodes or satellite passages triggering star formation can be done in future work.

Current-day element abundance measurements are reported with uncertainties of ≈0.02–0.05 dex (Ahumada et al. 2020; Jönsson et al. 2020). Consequently, we examine how the clustering changes once we incorporate errors in the chemical abundances and how this impacts our ability to trace birth properties with observational data.

For each star in both Chem2d and Chem15d , we redraw a new set of element abundances, each from a Gaussian distribution where the standard deviation is representative of the measurement uncertainty. We test two precision regimes: σ_err = 0.02 and 0.05 dex. The left column of Figure 8 shows the overlap in the age–R_birth plane of groups defined in Chem2d (top) and Chem15d (bottom) when the abundances are modified with an equivalent of 0.02 dex error in each abundance direction (black points and line). With the current best observational error of σ_err = 0.02 dex, the groups defined in both the two-dimensional and 15-dimensional abundance space retain separate birth properties, similar to the overlap found when the simulations have no error added (gray points and line). When abundances are redrawn with observational errors of σ_err = 0.05 dex for each data point (right column), we find that the majority of groups from the higher-dimensional simulation retain more separate birth information compared to those found in Chem2d.

Zoom In Zoom Out Reset image size

While current large surveys have captured many millions of stars, we have so far examined only about 30,000 stars with precise abundance measurements, within a narrow evolutionary state (e.g., apogee DR14 RC catalog; Bovy et al. 2014). To test more generally how useful chemical abundances are for linking to birth properties, we need to examine the impact of sample size.

So far in this work, we have been working with ∼229,000 and ∼44,000 star particles for the Chem2d and Chem15d simulations, respectively. Now we examine how the groups defined in abundance space, both with and without the addition of errors, change in the age–R_birth plane for a random subsample of 30,000 stars throughout the whole disk. Figure 9 shows the percent overlap for both the subsampled Chem2d (red) and subsampled Chem15d (blue) with no error (left), the equivalent of 0.02 dex error (middle) and 0.05 dex error (right) in each dimension. In order to test the consistency of these results, we subsample with 50 replications. The mean and standard deviation of the median percent overlap are shown as a point and ribbon.

Zoom In Zoom Out Reset image size

We can see from Figure 9 that as error increases, the amount the Chem2d groups overlap in the age–R_birth plane also increases. This indicates that for the sample size and error used in paper I of this series (Ratcliffe et al. 2020), only two abundance dimensions are not enough to recover separate birth time and place groups observationally for more than ≈6 groups.

On the other hand, the groups in the subsampled Chem15d are less affected by error than the subsampled Chem2d . We can see that as error increases, the amount each group overlaps with other groups in the age–R_birth plane stays more consistent than the groups in the subsampled Chem2d. While the groups from both simulations trace less birth information in the presence of subsampling, we see that the recovery of separate birth place and time groups is still possible with the inclusion of more abundances and nucleosynthetic families, especially for a smaller number (≤≈10) of groups.

The simulations used in this work (g7.55e11 with two abundances and g8.26e11 with 15 abundances) are Milky Way analogs from the NIHAO-UHD suite. Both simulations were bulge-dominated systems up to redshift z ≥ 1 with prominent stellar disks forming 7–8 Gyr ago. The formation of the α sequences was due to a gas-rich merger, with the high-α sequence forming during the early galaxy and the low-α sequence forming after the merger. The main differences between these simulations are a slightly different formation history (sampling valid formation histories of Milky Way–like galaxies) as well as an updated chemical enrichment prescription for the Chem15d model galaxy. The modifications made to chemical enrichment prescriptions are described in Buck et al. (2021) and enabled us to follow 15 different elements while at the same time leaving global galaxy properties such as star formation history, stellar mass, and disk size unaffected. We believe that these simulations are representative of the Milky Way due to their formation history and age gradient in the abundance plane.

Even though simulations provide us with particles and not individual stars, they allow us to examine the relationship between observable chemical properties, age, and birth location. Thus we focus on disk particles in our work.

In this work, we focused on three ways of grouping stars: hierarchical clustering, EnLink, and binning in the [α/Fe]–[Fe/H] plane. Section 4.4 shows that binning in just ([Fe/H], [α/Fe]) does not effectively link abundance information to birth properties, especially when there is higher-dimensional abundance data available.

Of the clustering methods we explore, hierarchical clustering is advantageous for observational work. Section 4.1.1 shows that leveraging density with an adaptive distance metric in Chem2d is the best way for chemical groups to correspond to distinct groups in the age–R_birth plane. However, due to the curse of dimensionality, leveraging density in a 15-dimensional space is difficult and unrealistic. We attempted to run EnLink in the 15-dimensional chemical abundance space of Chem15d; however, the clustering results were inconclusive and did not define many stars to a cluster. Additionally, for EnLink (or any other density-based clustering method) to be used correctly on Milky Way catalogs, the survey-selection function would need to be accounted for, as the selection function would possibly alter the distribution of stars in abundance space. Furthermore, we found that EnLink performed poorly when the density structure in the two-dimensional abundance plane of Chem2d vanished after the addition of observational uncertainty.

In this paper, we have focused on the relationship between abundances and birth properties of stars when the chemical bimodality is caused by a merger and successive dilution of the interstellar medium. We also consider the relationship when a galaxy is formed by clumpy star formation ⁸ . We examine the N-body and smooth-particle hydrodynamics simulation of the formation of an isolated galaxy outlined in Beraldo e Silva et al. (2021). Star-forming clumps at high redshift start forming low-α stars, then quickly self-enrich in α-elements due to their high-star-formation-rate density and produce a high-α sequence while a low-α sequence is produced by radially distributed star formation. After about 4 Gyr, the clumps become less efficient, and the high-α sequence stops growing (Clarke et al. 2019). For more detailed information of the simulation, see Beraldo e Silva et al. (2021) and Fiteni et al. (2021).

Figure 14 shows the [α/Fe]–[Fe/H] plane for a 230,000 particle subsample of simulation M2_c_nb, which undergoes clumpy star formation. Similar to Chem2d and Chem15d, there is a linear trend between the abundances and R_birth, where R_birth decreases as [Fe/H] increases. Age, however, does not appear to have a simple relationship between the two abundances. For instance, age decreases as [α/Fe] decreases for solar metallicity stars at higher values of [α/Fe], but the relationship is reversed for lower [α/Fe].

Zoom In Zoom Out Reset image size

We find that in this simulation, age and [Fe/H] are able to predict R_birth within ±0.72 kpc, about 40% better than the precision for the simulations focused on in this work. Again, we also see that the addition of other abundances does not notably change our ability to estimate R_birth, where the accuracy only increases by 0.01 kpc when [α/Fe] is included in the regression. This shows that the formation history in all three simulations sets an underlying relationship with age, [Fe/H], and R_birth, where if we know the metallicity and age of a star, we can determine where it was born.

However, while [Fe/H] and age are a link to R_birth, in this particular simulation the star formation history gives rise to a more complex relationship between the abundances ([Fe/H], [α/Fe]) and age (Figure 14). Therefore, chemically similar groups of stars identified using hierarchical clustering no longer correspond to separate groups in the age–R_birth plane in this scenario.

In order to determine the ability to extend our conclusions to the Milky Way, we compare the Milky Way's age and age dispersion in the [α/Fe]–[Fe/H] plane to the three simulations with simulated "observational" ages by redrawing from a normal distribution with their true age as the mean and a standard deviation of ≈2.6 Gyr, the median uncertainty of low-α stars from the Lu et al. (2021) catalog. When simulating observational ages for Chem2d and Chem15d , the dispersion in age across [α/Fe]–[Fe/H] is uniformly low (see middle rows of Figure 15), while M2_c_nb has an increase in age dispersion as [α/Fe] decreases (see bottom row of Figure 15). The Milky Way (shown in the top row of the same figure using ages and abundances from Lu et al. 2021) has a consistently small dispersion in age of about 2–3 Gyr, with the dispersion being slightly smaller for low ([Fe/H], [α/Fe]), the reverse of M2_c_nb's dispersion trends. On the other hand, similar to the Chem2d and Chem15d simulations focused on in this paper, the age of the stars in the Milky Way increases as [α/Fe] increases for a given value of [Fe/H].

Zoom In Zoom Out Reset image size

The question of which simulation most closely matches the Milky Way's star formation history still requires further investigation. The selection cuts described in Section 2.1 produce different density trends in the [α/Fe]–[Fe/H] plane for the three simulations and observational data, where some samples have both high- and low-α stars (e.g., Chem15d ), while others primarily consist of stars with lower values of [α/Fe] (e.g., M2_c_nb). There is room for exploration into how the Milky Way compares to the different [α/Fe]–[Fe/H] trends each simulation produces with different selection cuts; however, we find that our results are consistent under different cuts to capture disk stars.

This exploration shows that resolving distinct birth properties of chemically similar stars requires a small dispersion in age and R_birth trends in the [α/Fe]–[Fe/H] plane. Since the Milky Way has been shown to have age trends in the [α/Fe]–[Fe/H] abundance plane with small dispersion, we believe the conclusions we have drawn in this paper are relevant to the Milky Way. We conclude that the six groups found in our previous work (Ratcliffe et al. 2020) using hierarchical clustering in the 19-dimensional apogee red clump sample are expected to have (R_birth, age) distributions that differ from each other and reflect a continuous evolution of the disk.

There are some limitations to our approach. In Section 5.2 we discuss the effect of subsampling data with observational errors. However, we do not take into account the complexity of survey-selection functions. Additionally, in Section 5 we explore how errors affect our clustering results by adding 0.02 dex or 0.05 dex error to every abundance. Realistically, though, some abundances are measured more accurately than others.

While in the previous section we argued that the main conclusions of this work could be extended to the Milky Way and that the groups found in Ratcliffe et al. (2020) are expected to occupy different regions in (age, R_birth), we do not extend the numerical relationship between birth properties and abundances (Section 6) to estimate age and R_birth for stars in the Milky Way. This extension would require the assumption that the regression coefficients of the Milky Way are exactly the same as those used for the simulations. However, abundance values between simulations and the Milky Way differ as well as spatial coverage and orbital and structural properties. Our analysis in Section 7.3 shows that a quadratic regression model will probably work for the Milky Way, but determining the coefficients that are appropriate for the Milky Way is beyond the scope of this work.

Section 7.2 discussed how density-based clustering failed when a structure vanished due to measurement error or in a high-dimensional abundance space. For future work, combining hierarchical clustering with an adaptive distance metric would be interesting to explore. Partnering an adaptive distance metric with hierarchical clustering would avoid the problems of estimating density and determining the best distance metric in a high-dimensional-abundance space, and could potentially provide even more striking results.

As simulation resolution continues to increase, in the future it would also be useful to explore if satellite debris could be picked up by abundance clustering and complement clustering analysis done in action space (such as in Wu et al. 2022). In the simulations we use in this study, only ≈20 stellar particles have R_birth ≥ 20 kpc. With future data sets and simulations in mind, testing to see if accreted material differ in abundance space could be useful to determine accreted debris in the Milky Way.