Open Cluster Chemical Homogeneity throughout the Milky Way

Earlier studies have highlighted the possibility that the uncertainties for some of the APOGEE abundances in earlier data releases are overestimated (e.g., Ness et al. 2018), and our own preliminary analyses of the abundance dispersions in our clusters supported this assessment. Prior to DR16, uncertainties on the abundances in the APOGEE data releases were determined by examining the spread of abundances in well-sampled clusters (Holtzman et al. 2018), which required assuming that the cluster had no intrinsic spread. This assumption makes the uncertainties unsuitable for the study of chemical homogeneity.

To address this, we recalculate the random abundance uncertainties for the stars in our analysis using an improved method. ¹³ We worked closely with the APOGEE team to subsequently adapt our approach into the DR16 uncertainty determination. The approach used in DR16 (Jönsson et al. 2020) relies on a parametric fit to calculate the uncertainties for all the APOGEE stars. In this paper, however, we adopt a nonparametric approach because we find that a simple analytic function cannot adequately capture the relationships between uncertainty and signal-to-noise ratio (S/N), T_eff, and [M/H].

In brief, we use the differences in [X/Fe] values derived by ASPCAP for independent spectra of the same star to compute a relationship between the abundance uncertainties and S/N, T_eff, and [M/H] (Section 2.2.2). Then, we use this relationship to compute the uncertainties for our cluster members (Section 2.2.3).

In general, the ASPCAP pipeline is run on the stacked spectrum of each star, which comprises all visits to that star. For a subset of APOGEE stars, however, ASPCAP is run on the individual visit-level spectra, providing multiple independent sets of stellar parameter and abundance measurements for single stars. ¹⁴ We use these sets to estimate the random uncertainty of the ASPCAP measurements as a function of S/N, ${T}_{\mathrm{eff}}$, and [M/H]. The variations between ASPCAP values for spectra of the same star, at the same S/N, provide a more realistic representation of the random measurement uncertainties than the ones derived from cluster dispersions in earlier APOGEE DRs.

We define a sample of stars, hereafter called the Uncertainty Training (UT) sample, with ASPCAP solutions derived from two or more visit spectra with similar S/Ns (Δ(S/N)/(S/N) ≤20%), resulting in similar temperature (${\rm{\Delta }}{T}_{\mathrm{eff}}\leqslant 100$ K) and metallicity (${\rm{\Delta }}[{\rm{M}}/{\rm{H}}]\leqslant 0.07$ dex) values. These similarity criteria are imposed to ensure that differences in [X/Fe] are not due to different global spectral fits. We restrict our analysis to giant stars (using the ASPCAP_CLASS column and a limit of log g < 3).

The final UT sample of 8729 stars is then divided into five bins of S/N: 50–70, 70–100, 100–130, 130–200, and >200. These bins are chosen to provide finer sampling at lower S/N, where the effect of the S/N on output parameters is larger; above ${\rm{S}}/{\rm{N}}\approx 150$, there is a weak relationship with uncertainties. We explored dividing the sample into finer bins in S/N. However, there were no significant changes in the final derived uncertainties, and there was an increased risk of undersampling each bin so that the estimated uncertainty would not capture the effect of the varying range of stellar parameters.

We adopt a Voronoi binning scheme in T_eff–[M/H], within these fixed S/N ranges, to ensure both reliable measurements of the uncertainties in less populated regions of the parameter space and high-resolution measurements where possible. The Python package vorbin (Cappellari & Copin 2003) is used to group the UT sample into 2D bins of T_eff and [M/H], targeting at least 30 stars per bin. The final bins are populated with between 9 and 69 stars per bin, with an average of around 33.

The differences between pairs of visit-level [X/Fe] values for individual stars can be used to compute the standard deviation of the distribution from which the pairs were originally drawn. We assume this distribution to represent the intrinsic ASPCAP random uncertainty. The quantities are related by

$$\begin{eqnarray}&&{{\rm{e}}}_{[{\rm{X}}/{\rm{Fe}}],{\rm{k}}}=\displaystyle \frac{\sqrt{\pi }}{2}\mathrm{median}(| {\left[{\rm{X}}/\mathrm{Fe}\right]}_{i}-{\left[{\rm{X}}/\mathrm{Fe}\right]}_{j}| ),\end{eqnarray} \tag{ 1 }$$

where ${e}_{[{\rm{X}}/\mathrm{Fe}],{\rm{k}}}$ is the abundance uncertainty associated with the kth bin of T_eff and [M/H] at a given S/N, and ${[{\rm{X}}/\mathrm{Fe}]}_{i}$ and ${[{\rm{X}}/\mathrm{Fe}]}_{j}$ refer to abundance measurements derived from two independent visit spectra of the same star. The top row of Figure 1 shows the distribution of UT stars in [M/H] and T_eff, in bins of S/N, colored by the ${e}_{[\mathrm{Mg}/\mathrm{Fe}],{\rm{k}}}$ values. This demonstrates the complex pattern of ${e}_{[\mathrm{Mg}/\mathrm{Fe}],{\rm{k}}}$ derived in this way and the difficulty of describing such behavior with simple analytical expressions. Similar figures for all elements we analyze in Section 5 can be found in Appendix A (Figures A1 and A2).

Zoom In Zoom Out Reset image size

We looked for potential dependencies of ${e}_{[{\rm{X}}/\mathrm{Fe}],{\rm{k}}}$ on [X/Fe] itself—e.g., that stars with enhanced [Mg/Fe] may have different uncertainties than stars with solar [Mg/Fe] at the same ${T}_{\mathrm{eff}}$, [M/H], and S/N—by deriving ${e}_{[{\rm{X}}/\mathrm{Fe}],{\rm{k}}}$ in the same manner as above, but separately for stars with high and low [Mg/Fe]. We found no significant differences in the ${e}_{[{\rm{X}}/\mathrm{Fe}],{\rm{k}}}$(${T}_{\mathrm{eff}}$, [M/H], S/N) patterns. We also found no difference in the results discussed in Section 5 when using uncertainties derived from the solar [Mg/Fe] part of the UT sample compared to when using the full sample.

Given the computed array of ${e}_{[{\rm{X}}/{\rm{Fe}}],k}$ for each bin of T_eff, [M/H], and S/N, we can then sort any other star into a bin and assign it ${e}_{[{\rm{X}}/{\rm{Fe}}],k}$ values. We perform this sorting by training a Gaussian naive Bayes classifier algorithm (Pedregosa et al. 2011) on the UT sample and then applying the trained classifier to our cluster member sample (Section 3.3). This results in a bin membership probability for each star; we assign each star an uncertainty summed from all of the Voronoi bins and weighted by each bin's probability for that star:

$$\begin{eqnarray}&&{e}_{[{\rm{X}}/{\rm{Fe}}]}=\displaystyle \sum _{k}{p}_{k}{e}_{[{\rm{X}}/{\rm{Fe}}],{\rm{k}}},\end{eqnarray} \tag{ 2 }$$

where k is the bin index, p_k is the probability of bin k for this star, and ${e}_{[{\rm{X}}/\mathrm{Fe}]}$ is the final uncertainty for the star. Adopting this weighted average uncertainty ensures that all stars falling within a range of stellar parameters will not be assigned identical values for their abundance uncertainties due solely to the binning scheme. This approach also smooths the transitions between the bin edges. Ultimately, the uncertainties assigned to a star with stellar parameters near a bin edge are not driven by the distribution of the bins but by the uncertainties estimated from the UT sample that fall within a neighboring area of the star in [M/H] and T_eff.

The bottom three rows of Figure 1 show examples of these rederived uncertainties for [Mg/Fe], [Ni/Fe], and [Na/Fe] (note the differences in color scaling) along [M/H] and T_eff, in bins of S/N. All show the expected improvement in precision with higher S/N. Other generic patterns are also clearly visible—for example, the increase in uncertainty at low metallicities and/or high temperatures, where lines become weaker, and at very low temperatures, where lines become increasingly blended. Mg and Ni have lower uncertainties compared to Na, which is expected due to the difficulty in measuring Na lines in APOGEE. Each element also has its own unique patterns, reflecting the range of difficulty in measuring lines of different elements in different parts of stellar parameter space.

The DR16 uncertainties are not systematically higher or lower than the uncertainties derived here for stars in our cluster sample (Section 3), and we find qualitatively similar results for the analysis in Section 5 if we use either set of values.

For each cluster, we start our membership search with APOGEE and Gaia-DR2 stars within twice the cluster radius, 2R_cluster, using central coordinates and "total" cluster radii (their r₂ = R_cluster) from K13 (Section 2.5). The "central stars" (within 1R_cluster) define the cluster's kinematical signature, and the "annulus stars" (between 1.5–2R_cluster) define the background distribution (e.g., Figures 3(a)–(b)). We only consider the 366 K13 clusters that have six or more stars within 1 R_cluster in the APOGEE catalog that meet the quality criteria above.

We search for RV peaks associated with comoving stars in each cluster (e.g., Figure 3(c)) location by subtracting a kernel density estimate (KDE) of the annulus stellar RV distribution (shown in green) from that of the central stars (shown in blue). The residual (shown in red) is then fit with a Gaussian to determine the central RV ($\langle \mathrm{RV}\rangle$) and the width (σ_RV) of the dominant peak. We also measure the ratio of the Gaussian amplitude (A_RV) to the standard deviation of the residuals (σ_resid) more than 3σ_RV away from the Gaussian center; this metric quantifies the strength of the RV signal. Visual inspection demonstrates that for our clusters, an A_RV/σ_resid > 9.5 corresponds to a cluster in RV space. Smaller values tend to be dominated by noisy residuals driven by a low number of annulus stars.

We obtain PM information for all stars within 2R_cluster using the Gaia TAP+ query from the astroquery package in Python and keep stars that pass the quality cuts mentioned in Section 2.3.

Next, we ensure that the Gaia-DR2 data sample has the same color–magnitude range as the APOGEE stars by matching to the "apogeeObject" files used in the APOGEE targeting pipeline, which contains the 2MASS (Skrutskie et al. 2006), Spitzer–IRAC GLIMPSE (Benjamin et al. 2005; Churchwell et al. 2009), and AllWISE (Wright et al. 2010; Cutri et al. 2013) photometry used to calculate extinction (Majewski et al. 2011; Zasowski et al. 2013b, 2017). We then restrict the Gaia-DR2 stars to the same ${(J-{K}_{s})}_{0}$ and H limits sampled by the APOGEE stars in the vicinity of that cluster (generally ${(J-{K}_{s})}_{0}\geqslant 0.5$ and $7\leqslant H\lt 12.2$).

We use this cut to ensure that the PM distribution obtained from Gaia-DR2 stars is an accurate representation of the APOGEE stars that we are considering for membership. However, we compare the membership with and without using this cut for all 10 clusters studied in Section 5. Although the PM distribution is altered slightly, the final cluster members determined are the same for these clusters irrespective of the color–magnitude cut.

In a small fraction of cases (6%), the apogeeObject files do not span the full background annulus region, but we have confirmed that the distributions of μ_R.A. and μ_decl. do not change across the small angular scales of our clusters, so we consider even these partial annuli to be representative of the background distributions.

As with the RVs (Section 3.1.2), we compare KDEs of the central and annular distributions to identify any signal of the cluster, this time with 2D KDEs (μ_R.A. × μ_decl, shown in Figure 3(d)). Because the annular PM distribution is much less noisy than in the case of RVs, we model the entire central PM distribution (shown in blue) as the sum of a 2D Gaussian and a scaled copy of the annular PM distribution (shown in green). The best-fit Gaussian (shown in red) center ($\langle \mathrm{PM}\rangle$) and 2D dispersion (σ_PM(R.A., decl.)) are taken as the PM distribution of the comoving cluster stars.

We first compute cluster membership probabilities, based on RVs and PMs, for each star within 2 R_cluster of a cluster. These probabilities are the values of Gaussian distributions with the means and standard deviations derived from the RV and PM fitting in Sections 3.1.2–3.1.3, scaled to have a maximum value of unity. We consider as likely members stars that fall within a 3σ window (shown in purple in Figures 3–5) on the combined RV–PM probability. We give equal weighting to each kinematic dimension while calculating the combined probability, because weighting each dimension on how distinct it is from the background did not yield any changes in the final selected cluster members. We choose a selection window of 2σ on the combined probability for the analysis in Section 5 because we observed a few outliers in the metallicity distribution (e.g., Figure 3(f)) that were removed when we used a stricter cut.

Examples of this entire procedure and its results are demonstrated in Figures 3–5. The color coding for all panels is as follows: locations and kinematical information for the central APOGEE stars are plotted in blue, for the annulus APOGEE stars in green, for the final members in purple, and for the Gaia-DR2 stars used for the PM background in gray points. The top row (panels (a) and (b)) show the on-sky distribution of stars, with the inner R_cluster in a blue circle and the outer annulus enclosed in green circles at 1.5 and 2 R_cluster. Panel (a) shows the stars used in characterizing the cluster (Sections 3.1.2–3.1.3), and panel (b) shows the final cluster members (Section 3.1.4).

The middle row (panels (c) and (d)) show the RV and PM distributions of the central and annulus stars, with the fitted residual peaks in red. The bottom row (panels (e) and (f)) are not used for membership selection and are only examined when setting reliability flags for entire clusters (Section 3.2).

Figure 3 shows the recovery of the well-studied cluster NGC 6819 (e.g., Hole et al. 2009; Platais et al. 2013; Yang et al. 2013; Wu et al. 2014; Lee-Brown et al. 2015), Figure 4 shows the recovery of the less-studied cluster FSR 0494 (Froebrich et al. 2007; Zasowski et al. 2013a; Donor et al. 2018), and Figure 5 shows the nonrecovery of the cluster cataloged as ASCC 166 (e.g., Kharchenko et al. 2005; Cantat-Gaudin et al. 2018).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We calculate the 3D velocity dispersion of a cluster, a proxy for cluster mass, from its RV and PM dispersions and heliocentric distance using the following equation:

$$\begin{eqnarray}&&{\sigma }_{{\rm{tot}}}=\sqrt{{\sigma }_{{RV}}^{2}+({\sigma }_{{\mu }_{\alpha }}^{2}+{\sigma }_{{\mu }_{\delta }}^{2}){d}_{{\rm{helio}}}^{2}},\end{eqnarray} \tag{ 4 }$$

where σ_tot is the space velocity dispersion, corrected by the uncertainties as described in Section 4; ${\sigma }_{{RV}},{\sigma }_{{\mu }_{\alpha }},$ and ${\sigma }_{{\mu }_{\delta }}$ are the dispersions in the cluster for each kinematic dimension; and d_helio is the heliocentric distance, assigned as the median of the stellar member distances from RA17 (Section 2.4). We observe a strong correlation between the calculated RV and PM velocity dispersions, which ensures that σ_tot is not being driven by any one dimension alone.

We find that the cluster chemical scatter is positively correlated with the space velocity dispersion of the cluster at relatively high levels of significance (p < 0.019) for [Mg/Fe] and [Ca/Fe]; possibly positively correlated at a low level of significance (0.038 < p < 0.059) for [Ni/Fe], [Si/Fe], [Al/Fe], and [Fe/H]; and uncorrelated (p > 0.38) for [Cr/Fe] and [Mn/Fe]. Figure 9 shows the intrinsic scatter in [Fe/H] and the rest of the abundances as a function of space velocity dispersion, along with the associated Spearman correlation coefficients (C) and p values.

Zoom In Zoom Out Reset image size

To understand why only certain elements show this trend between intrinsic scatter and σ_tot, we look for natural ways to group elements based on their properties. For example, we notice that this trend is not exclusive to the α elements that we study. Although intrinsic scatter is positively correlated with σ_tot in Mg, Ca, and Si (albeit at low significance), we observe a similar trend in an odd-Z element (Al) and an iron-peak element (Ni) at lower significance. In Section 5.3 we discuss how the abundance scatters we observe behave for elements that show abundance variations and anticorrelations in GCs (e.g., Mg, Al, and Si). Again, this trend is not restricted to these three elements. However, we find that if we group elements based on their dominant nucleosynthetic process of production, we see a distinction between elements that are produced predominantly by core-collapse supernovae (CCSNe) versus Type Ia SNe.

In order to visualize the differences in the strengths of correlation in groups of elements and to explore the roles that different enrichment events may have played in them, we use the empirically determined fractional contribution of CCSNe (f_cc ) for each element from Weinberg et al. (2019). Figure 10 shows the Spearman correlation coefficients from Figure 9 against f_cc , with the points colored by the p value of their correlation. Here, f_cc represents the fraction of each element contributed by CCSNe at a given metallicity [Mg/H], assuming that these elements are produced exclusively by Type Ia SNe and CCSNe. We calculate f_cc from Equation (11) in Weinberg et al. (2019), using the median [Mg/H] value for each of our clusters.

Zoom In Zoom Out Reset image size

From Figure 10, in addition to metallicity (traced by [Fe/H]), the elements that show a correlation between intrinsic abundance scatter and cluster velocity dispersion (with C > 0.6) at both higher (with p < 0.019) and lower (0.038 < p < 0.059) levels of significance are the ones that are produced mostly by CCSNe.

While these correlations between abundance scatter and space velocity dispersion are interesting, we note that the statistical analysis is done using only 10 clusters and the elements that we list as correlated have varying levels of significance (p values). Here we explore the caveats associated with these correlations. One effect of the low sample size is seen when we randomly remove any one cluster from the analysis using a jackknife resampling. For certain elements, the correlation becomes insignificant if we remove a specific cluster (e.g., removing NGC 2204 in the case of Si or removing NGC 6791 in the case of Ni) from the analysis. However, we found that no one cluster is responsible for systematically reducing the significance across the set of all the elements. Because we arrive at this cluster sample by preferentially selecting clusters with high quality membership and reliable abundance uncertainties (Section 5.1), we have no a priori reason to drop any particular cluster from any particular elemental trend.

We explored the effect of using a different correlation metric (Kendall's tau) on the strength and significance of the correlations seen in Figure 9. We do not use the Pearson metric in this comparison because it is susceptible to outliers and assumes linear relationships. The Spearman and Kendall correlation metrics are highly correlated themselves, though their magnitude is not equal, and they are more robust to outliers as they use the ranks of the variables rather than the actual values. We show the Spearman metric in Figures 9 and 10 because it is more commonly used, although we note that Kendall may be more accurate for small sample sizes. The significance of the correlations for Al, Si, Fe, and Ni depends somewhat on the choice of the metric, with Kendall's tau associated with p values up to 0.07, which reinforces our classification of these elements as possible correlations. However, Mg and Ca scatters are significantly correlated with velocity dispersion irrespective of the metric used, while Mn and Cr remain uncorrelated.

As discussed in Section 5.5.1, it is interesting that the elements that show some level of correlation with space velocity dispersion (Figure 9) are those that are predominantly produced by CCSNe. However, due to the small sample size, both in clusters and elements, and other caveats discussed in Section 5.5.2, we cannot definitively conclude that these correlations are significant. However, in case these findings are validated by future larger studies, here we explore the potential implications of this result.

Figure 10 suggests that the nucleosynthetic processes that are responsible for the production of elements in the ISM may have an observable effect on the final abundance scatter within the cluster. It also hints at the existence of a difference in the ejection radii between the pollution mechanisms of Type Ia SNe and CCSNe.

Why does the chemical scatter within a cluster depend on the mass of the cluster? Equally interesting, why is this correlation present in certain elements and not in others? The dependence of the intrinsic abundance scatter on cluster mass can be understood by looking at what we know about cluster formation processes. Clusters are formed from giant molecular clouds in filamentary structures (Kounkel & Covey 2019), forming strings of star-forming gas. So, gas accreted to form more-massive clusters will not be limited to a sphere surrounding the final cluster, but rather spans a larger range in distance. Fujii (2015) explored the possibility that young massive clusters may be formed by hierarchical merging of subclusters or smaller open clusters.

Although the OCs we are using in our analysis are not as massive as the young massive clusters discussed in Fujii (2015), this mechanism hints that more-massive clusters that we study may have been formed by accreting gas over larger ranges in distance in the initial cloud. This would result in the more-massive clusters having a larger scatter in metallicity and abundances of certain elements depending on their correlation lengths and mixing efficiency in the initial clouds before star formation began. Krumholz et al. (2019) also suggest that massive clusters are formed over extended formation times rather than a single freefall timescale. This could again increase the chances of the star-forming cloud being polluted by exploding high-mass stars that have already formed.

How are variations in the correlation lengths of different metals in an initial cloud reflected in the chemical homogeneity of the final clusters formed from the cloud? Armillotta et al. (2018) studied this using hydrodynamical simulations, and they observed that if the initial field of the metal is correlated over smaller distances (<6 pc), turbulent mixing will efficiently smooth out these inhomogeneities. Metal fields that are correlated over much larger distances (>40 pc) will also yield homogeneous stellar chemistry as the variations span a range larger than the typical cloud size. However, metals that are correlated on intermediate lengths in the ISM (6–40 pc) can have a higher scatter in the stellar abundances of their final clusters. So, with this reasoning, elements that have correlation lengths within the intermediate range in the initial cloud may be expected to have a higher abundance scatter in the final stellar members for massive clusters.

Furthermore, because elements belonging to different nucleosynthetic groups have been shown by Krumholz & Ting (2018) to have different correlation lengths in the initial cloud, we may observe this trend only in certain elements or nucleosynthetic groups. They propose that there should be no significant differences between the correlation lengths of Type Ia SNe and CCSNe since both types of explosions have comparable energy budgets. However, Figure 10 suggests the presence of a quantitative difference between the correlation lengths of these two mechanisms that pollute the ISM, and that this difference may manifest itself in the abundance scatter of nucleosynthetic element groups in the final stellar populations.