A Comprehensive Study of Open Cluster Chemical Homogeneity using APOGEE and Milky Way Mapper Abundances
Amaya Sinha Gail Zasowski Peter Frinchaboy Katia Cunha Diogo Souto Jamie Tayar Keivan Stassun


0009-0005-0182-7186]Amaya Sinha
Department of Physics & Astronomy,  University of Utah, 115 South 1400 East, Salt Lake City, UT 84112, USA

0000-0001-6761-9359]Gail Zasowski
Department of Physics & Astronomy,  University of Utah, 115 South 1400 East, Salt Lake City, UT 84112, USA

0000-0002-0740-8346]Peter Frinchaboy
Department of Physics & Astronomy, Texas Christian University, Fort Worth, TX 76129, USA

0000-0001-6476-0576]Katia Cunha
Steward Observatory, University of Arizona, Tucson, AZ 85721, USA
Observatório Nacional, Rua General José Cristino, 77, 20921-400 São Cristóvão, Rio de Janeiro, RJ, Brazil

0000-0002-7883-5425]Diogo Souto
Departamento de Física, Universidade Federal de Sergipe, Av. Marcelo Deda Chagas, Cep 49.107-230, São Cristóvão, SE, Brazil

0000-0002-4818-7885]Jamie Tayar
Department of Astronomy, University of Florida, Bryant Space Science Center, Stadium Road, Gainesville, FL 32611, USA

0000-0002-3481-9052]Keivan Stassun
Department of Physics and Astronomy, Vanderbilt University, VU Station 1807, Nashville, TN 37235, USA

\begin@twocolumnfalse
    § ABSTRACT 
Stars in an open cluster are assumed to have formed from a broadly homogeneous distribution of gas, implying that they should be chemically homogeneous. Quantifying the level to which open clusters are chemically homogeneous can therefore tell us about ISM pollution and gas-mixing in progenitor molecular clouds. Using SDSS-V Milky Way Mapper and SDSS-IV APOGEE DR17 abundances, we test this assumption by quantifying intrinsic chemical scatter in up to 20 different chemical abundances across 26 Milky Way open clusters. 
We find that we can place 3σ upper limits on open cluster homogeneity within 0.02 dex or less in the majority of elements, while for neutron capture elements, as well as those elements having weak lines, we place limits on their homogeneity within 0.2 dex. Finally, we find that giant stars in open clusters are ∼0.01 dex more homogeneous than a matched sample of field stars.
    
\end@twocolumnfalse



§ INTRODUCTION


Immediately after the Big Bang, the only elements in the universe were hydrogen, helium, and trace lithium. It took the formation of stars and galaxies to populate the universe with the rest of the periodic table. Therefore understanding where and how stars produce and disperse heavy elements is essential to understanding the enrichment of the universe. However, many questions regarding the chemical enrichment of the universe still remain unanswered. Specifically, there is still much uncertainty on how well-mixed giant molecular clouds are or how heavy elements get from their production sites into stars <cit.>. Fortunately, with a few exceptions, the surface abundances of stars are fossil records of the gas composition from the molecular cloud in which they formed. As a result, we can use the present-day chemistry of stars to learn about the chemistry of the Milky Way in the past. 

In the age of large astronomical surveys such as GALAH <cit.>, LAMOST <cit.>, RAVE <cit.>, APOGEE <cit.>, and Gaia <cit.>, we can now probe the chemistry of stars in the Milky Way on the scale of ∼ 0.1 dex or smaller in multiple elements across different nucleosynthetic families, allowing us to trace different chemical enrichment pathways. Furthermore, we can now study the chemistry of the Milky Way at multiple different scales, from the simplest population in conatal binaries <cit.> to large populations of dispersed field stars <cit.>. 

Stars in an open cluster (OC) are assumed to have formed from a broadly homogeneous distribution of gas at the same time, implying that they should all have the same age and be at the same distance <cit.>. Using the chemistry of OC stars, we can infer the chemistry of the gas available at that point in the Milky Way's history, in particular within the thin disk. In the past, it has been suggested that using assumptions of chemical homogeneity from simple stellar populations like open clusters, it would be possible to reconstruct a dissolved cluster purely by its members' chemistry <cit.>. This technique, known as chemical tagging, has been a strong motivator for studies of cluster chemistry.

While many studies support this assumption of OC chemical homogeneity <cit.>, there has been work showing that at least some clusters are chemically inhomogeneous <cit.>. For example, <cit.> argued that NGC 6791 may not be chemically homogeneous, due to the presence of a potential Na-O anti-correlation, a relationship most commonly found in globular clusters <cit.>. This would be an exciting result as NGC 6791 already unique in the Milky Way, as both the most massive and the most metal rich open cluster. However other studies, such as <cit.>, have shown that it is chemically homogeneous within measurement uncertainties. For a more detailed discussion on previous studies of open cluster chemical homogeneity, see Section <ref>.

There are reasons why OCs could be heterogeneous in specific elements. Slow neutron capture element abundances, such as Sr, Ba, and Zr, can change over a star's lifetime as it enters the AGB phase of its evolution <cit.>. Dredge-up, which occurs in stars on the giant branch, causes the star's convective envelope to expand, and it eventually gets deep enough to pull CNO-cycled elements to the surface, thereby altering the surface abundances we measure <cit.>. NGC 6705 is an interesting OC regarding this effect, as it has been observed to also be enhanced in Na due to dredge up <cit.>.

The surface abundances of elements such as Mg can be affected by effects like mass transfer and atomic diffusion <cit.>. However, the latter only weakly impacts the upper giant branch. Lastly, mass transfer <cit.> and pollution events such as planetary engulfment <cit.> 
can also alter a star's surface abundances.

However, if an open cluster was measured to have nonzero chemical scatter even accounting for these factors, that could point to interesting and understudied physics that may have occurred during the formation of the OC. Simulations have shown that turbulent mixing during cloud assembly naturally produces a stellar abundance scatter that is ∼ 5 times smaller than that in the natal gas <cit.>; suggest that this mixing could explain the observed chemical homogeneity of stars forming from the same molecular cloud. This is supported by recent work by <cit.> who find open clusters in FIRE-2 simulations <cit.> to have chemical scatter within 0.02 dex on average.

However, chemical inhomogeneity in real clouds could be due to effects not fully captured by simulations, related to internal turbulence and gas mixing within the progenitor molecular cloud or pollution events such as core collapse supernovae (CCSNe) that occurred earlier in the cluster's lifetime <cit.>. 


Quantifying the level of chemical homogeneity in open clusters across a broad set of elements from various nucleosynthetic families would provide the basis for understanding the physics of early OC formation. 

This work aims to constrain the chemical homogeneity in a large set of abundances and clusters to disentangle the causes of those chemical variances. The structure of the paper is as follows: Section <ref> outlines the survey data, verification of the abundance uncertainties, and determination of the cluster membership. Section <ref> details the methodology and calculation of the intrinsic scatter within each [X/Fe] across the final cluster sample. Section <ref> presents the results of our work, and Section <ref> compares our results to previous findings.




§ DATA





 §.§ SDSS




  §.§.§ SDSS-V/MWM

The abundances and radial velocities (RVs) we use are primarily drawn from the Milky Way Mapper (MWM; J.A. Johnson, in prep), a component of the fifth generation of the Sloan Digital Sky Survey <cit.>. We use data from Internal Product Launch 3 (IPL-3), which will form the basis for SDSS DR19 (K. Hawkins, in prep). This dataset builds off of
the observing strategies and survey goals outlined in SDSS Data Release 18 and includes observations of over a million targets <cit.>.

SDSS-V/MWM uses two telescopes: the Sloan Foundation Telescope at APO <cit.> and the duPont Telescope at LCO <cit.>. Both are outfitted with nearly identical custom-built 300-fiber APOGEE spectrographs <cit.>, which reach a resolution of R ∼ 22,500, spanning the range of wavelengths between 1.51-1.70 μm. Unlike in SDSS-IV, which used a plug-plate system, SDSS-V now uses robotic fiber positioners <cit.>, which benefited from the adoption of a three-element corrector for the Sloan telescope at APO <cit.>. 

Within IPL-3, three different data pipelines were used to analyze the data taken from APO and LCO: 
The Payne <cit.>, The Cannon <cit.>, and the APOGEE Stellar Parameters and Abundances Pipeline <cit.>. Both The Payne and The Cannon are label-transfer methods that determine stellar labels from spectra after being trained on a set of spectra with known labels. They are differentiated by the fact that The Cannon is a data-driven model which requires no information about stellar models. Rather measurements from the Cannon inherit information from the models from its training data. The Payne incorporates physical models directly into its analysis. While these datasets are similar in many aspects, comparing the limits derived through each of them will provide a stronger constraint on the true homogeneity of the OCs in our sample.



  §.§.§ SDSS-IV/APOGEE


In addition to IPL-3, we use abundances and RVs from the seventeenth and final data release <cit.> of SDSS-IV's <cit.> Apache Point Observatory Galaxy Evolution Experiment <cit.>, which contains over 700,000 stars. The initial targeting strategy for APOGEE-1 and APOGEE-2 are outlined in <cit.> and <cit.>, respectively, and the final targeting for APOGEE-2 is outlined in <cit.> for APOGEE-2 north, and <cit.> for APOGEE-2S. 

The details for the APOGEE data reduction pipeline are described in <cit.>, and the details for ASPCAP are found in <cit.>. The description for the updates to these pipelines for DR17 is included in Holtzman et al. (in prep). The MARCS model atmospheres and interpolation methodology used in APOGEE are described by <cit.> and <cit.>. The line lists used for DR17 are outlined in <cit.>, and the spectral fitting used for ASPCAP is described in <cit.>. The details describing the APOGEE spectral grids can be found in <cit.>, <cit.>, and <cit.>, and lastly, the details for Turbospectum can be found in <cit.> and <cit.>.

We also include abundances from the BACCHUS Analysis of Weak Lines in APOGEE Spectra <cit.>, a value-added catalog (VAC) in DR17. This VAC provided abundances for several chemical species having weak and blended lines that cannot be reliably analyzed using ASPCAP. This sample consists of high signal to noise (SNR > 150) red giant stars with no flags in either  or  and analyzed using the BACCHUS code <cit.>, which measures line-by-line elemental abundances from on-the-fly spectral synthesis. High quality measurements are stacked to create a sample of elemental abundances for elements with weak or blended lines. We use the BAWLAS VAC abundances and uncertainties for the following elements: Na, P, S, V, Cu, Ce, and Nd.

Two separate uncertainties are reported for each BAWLAS abundance measurement. One is the  describing the measured uncertainty from the combined spectra using the same methodology as ASPCAP. The other is , which describes the uncertainty derived from the spectral lines themselves. Here to remain consistent in our analysis we use  as it is the closest to the uncertainty calculation method that we verify in Section <ref>.

The calculation of the abundances in the BAWLAS catalog is outlined fully in <cit.>. Each spectrum has an associated , with values , , or . Both  and  indicate either suspicion with the final fit or total failure, and only  indicates that the spectral fit is trustworthy. To ensure the highest quality sample, we require all the stars used in this study to have measurements with  = 1.



  §.§.§ Quality Cuts


We limit our sample in APOGEE DR17 to stars with  0.1 km s^-1 and  5000 km s^-1 to ensure our stars have reliable radial velocities. We also restrict our sample to stars with  1 km s^-1, to remove potential binaries within our sample <cit.>. Here  refers to the average radial velocity derived from individual RVs that are drawn from cross-correlation of individual spectra with combined spectrum.  refers to the uncertainty on that radial velocity, and  refers to the scatter of individual visit RVs around the average. To ensure the sample has reliable measurements we enforce a . These limits are identical between DR17 and IPL-3. We also limit our sample to stars between 3000 K and 6500 K, and -1 ≤ [Fe/H] ≤ 1.

We also exclude from our sample any stars that have [X/Fe] flags in more than 2 elements. Lastly, we only sample from stars with  in order to ensure that every star we study is a member of the giant branch. These requirements result in a sample of 305,201 stars. When using IPL-3 we use the corresponding columns and limits, with the exception of  which is not included in IPL-3. 

In DR17's allStar file, we enforce quality cuts on the  and  columns, the details of which are included here: <https://data.sdss.org/datamodel/files/APOGEE_ASPCAP/APRED_VERS/ASPCAP_VERS/allStar.html>. The details of the APOGEE bitmasks are located here <https://www.sdss4.org/dr17/algorithms/bitmasks/>. Within the  column we enforced the following requirements:

 
  *  == 0; BAD overall for star: set if any of , , , , , SN error are set, or any parameter is near grid edge ( is set in any )
 
  *  == 0; FERRE failed to return a value for metals.
 
  *  == 0; Elemental abundance from window differs > 0.5 dex from parameter abundance for [α/Fe].

and within the  we made the following cuts:

 
  *  == 0; Star has neighbor more than 10 times brighter.
 
  *  == 0; Star has neighbor more than 100 times brighter.
 
  *  == 0; WARNING: cross-correlation peak with template significantly broader than autocorrelation of template.


The IPL-3 ASPCAP allStar file does not have a published  column. As such, we require that all the stars in our sample have  and . For all three IPL-3 pipelines we use we enforce a  requirement. Where possible we enforce a  requirement within each pipeline. Within the three IPL-3 allStar files we use, we enforce these quality cuts on the  column, which has bits that correspond to DR17's  column.



 §.§ Cluster Membership

We start from the catalog of cluster members published in <cit.>, hereafter , which contains membership information for over 200,000 stars across ∼ 2000 OCs using Gaia DR2 <cit.>. Due to the more recent availability of Gaia DR3 <cit.>, we first match the stars identified as cluster members in  to their DR3 kinematics and positions.  used two spatial (RA α, DEC β) and two kinematic (proper motion-RA δ _α *, proper motion-DEC δ _β *) parameters, as well as parallax ϖ, as inputs into an unsupervised machine learning algorithm to determine cluster membership. We limit the initial cluster membership candidacy to stars from  within three cluster radii of their cluster centers. Using each cluster's distribution in radial velocity, we find that a minimum probability cut of P ≥ 0.5 in the  catalog maximized the overlap in membership within DR17 and IPL-3 while also minimizing contamination from non-cluster members. Of the 2019 OCs identified in , only 145 clusters have any members within both DR17 and IPL-3.



  §.§.§ Kinematic Selection


To further ensure that the clusters identified had minimal contamination, we use a Kernel Density Estimator with a variable bandwidth, following Silverman's Rule, to measure the dispersion in four dimensions (radial velocity, δ _α *, δ _β *, and ϖ). We reject stars further than two standard deviations from the cluster median. For an example of this methodology, see Figure <ref>, which shows the final distributions after these cuts in M67. The kinematic selection plots for all clusters in our sample can be found in Appendix <ref>.

We plot the [Fe/H] distribution of each cluster against the distribution of nearby non-cluster members between 2–5 cluster radii, as published in CG2018. We visually inspect the HR diagrams for the cluster members compared to the annulus to ensure no contamination, as each cluster should follow a single distinct isochrone. We use MIST isochrones <cit.>, generated using ages from <cit.>, and the median cluster [Fe/H]. Lastly, we only select clusters with N ≥ 6 members, resulting in a final sample set of 26 open clusters. Determination of this minimum limit and the final sample size is related to material in Section <ref>.






  §.§.§ Final Cluster Sample

The distribution of the final cluster sample in radius, age, and [Mg/Fe] versus [Fe/H] can be seen in Figure <ref>. 

We calculate the positional, kinematic, and orbital information for each cluster in our final cluster sample. Using Astropy SkyCoords <cit.> we calculate the Cartesian X, Y, Z galactocentric coordinates for all the clusters in our sample, as well as the galactocentric radius.

We integrate each cluster's orbit to measure its Z_max, eccentricity, and guiding radius using Galpy, with the  gravitational potential <cit.> 
While this potential lacks a bar, no member of our sample is close enough to the Galactic center to produce a noticeable difference in the final measured eccentricity, guiding radius, or maximum height above the galactic plane. While there have been measured effects on these parameters due to the Milky Way's spiral arms, as seen in <cit.>, as we do not use these parameters in the final science results, we do not include spiral arms in our potential.

We use a Monte Carlo method with N = 100 iterations to estimate uncertainties on these parameters. We measure 3D space velocity dispersion as a proxy for cluster mass and using the methodology outlined in <cit.>, we quantify the ratio of nucleosynthetic enrichment within each cluster from CCSNe and Type Ia supernovae. 

Lastly, we include age estimates on all our clusters from <cit.>. These ages were derived using an artificial neural network trained on a sample of reference clusters. For further description see <cit.>. At this stage we also flag red clump stars by eye within each cluster. All the cluster parameters are included in Table <ref>, and will be included as a machine readable table.






 §.§ Verification of Abundances & Uncertainties


The values of intrinsic scatter in OC abundances are on the order of ∼0.01 dex <cit.>, as are the abundance uncertainties within DR17. Therefore, verifying that the uncertainties on the abundances we are studying were not underestimated or overestimated was necessary. The method we use is the same as <cit.>. Since each star has some intrinsic true abundance measurement, repeated observations of the same star should result in a normal distribution around that value, where the width of the distribution is due only to the measurement uncertainties. Therefore we can use stars that have multiple visits in different APOGEE fields to measure the true uncertainty, and compare it to reported values within DR17 and IPL-3. To quantify our multiple-visit empirical uncertainty, we use Equation <ref> from <cit.>.

    e_[X/Fe],k = √(π)/2median(|[X/Fe]_j-[X/Fe]_i]|),


where [X/Fe]_j and [X/Fe]_i are the abundance measurements from the same star's ith and jth visits. e_[X/Fe],k is the resulting [X/Fe] uncertainty after median stacking the pairwise measurements in the kth bin. 

We group stars by SNR in bins spanning 50–70, 70–100, 100–130, 130–200, and greater than 200, as shown in Fig <ref>. Within each SNR range, stars are binned by T_ eff and [M/H], where Δ[M/H] = 0.2 dex and Δ T_ eff = 200 K. We use the K-S statistical test to ensure the distribution of abundance differences in each bin is consistent with a normal distribution, and flag those that are not to ensure they do not contaminate our sample. Bins that are not consistent with a normal distribution have two main causes. Firstly, some are at the edge of the [M/H] parameter space, in particular at low metallicities where ASPCAP is less reliable. Secondly, some bins have poor completion, with less than 10 measurements.

Within those bins where the empirical uncertainties are well-measured and normally distributed, we find very good agreement with the native pipeline uncertainties; thus for the rest of this work, we adopt those pipeline uncertainties directly from DR17 and IPL-3. 



§ MEASURING CHEMICAL SCATTER




 §.§ Paired Stars Method


Our primary method of determining the intrinsic scatter uses the difference in abundances between stars close to one another on the HR diagram: Δ T_ eff < 100 K and Δlogg < 0.1 dex. These limits are chosen because the maximum induced abundance offset between pair members due to the systematics discussed in Section <ref> is on the order ∼0.001 dex.

We measure the intrinsic dispersion within each pair using Equation <ref>, where e_1 and e_2 are the abundance uncertainties of each star in the pair, |Δ[X/Fe]| is the absolute value of the difference in abundance measurements between the pair, and σ is the inferred intrinsic dispersion.


    σ = √(π/2|Δ[X/Fe]|^2 - (e_1^2 + e_2^2))/2)


We use a Monte Carlo method with N = 100 iterations to vary the abundance measurement of each star within the pair to estimate a final uncertainty on the measured pairwise scatter. Within each cluster, we separate the red clump and RGB stars as to not induce any scatter from potential evolutionary effects. Within each sub-sample, we sort the stars into pairs and then measure the intrinsic scatter between them. Finally, we take the median scatter of all the pairs within an OC as the true intrinsic scatter of the cluster. This method allows for extremely precise results, with final uncertainties on the order of ∼ 0.001 in most elements. 

To determine the number of pairs needed for a reliable measurement with this method, we use a synthetic “cluster” of points, with “true” and “observed” [X/H] abundances, and the same temperature and log(g) distribution as our real clusters. The true abundances for the synthetic stars reflect a given intrinsic scatter for the synthetic cluster. The observed abundances are generated by perturbing each true abundance by a random value drawn from a normal distribution with σ set to the uncertainty of a real APOGEE star with the same temperature, log(g), and metallicity.

We then pair the stars as described in Section 3.1 and perform the intrinsic scatter measurement described using N=3 to N=20 pairs. We find that the measured cluster scatters are noisy and have both systematic offsets and larger uncertainties up until N=8, at which point the difference between the true and measured scatter does not change at larger N. Thus, we require a minimum of eight stellar pairs for clusters using this method.

Due to the restrictions outlined above regarding the separation of the pairs in T_ eff and logg, as well as the minimum number of pairs required, this method can only be applied to 15 of the 26 clusters studied in this paper. However due to the significantly higher precision of these values as compared to those derived using the method outlined in Section <ref>, within these 15 clusters we only publish results from this method.




 §.§ Maximum Likelihood Estimator


We adopt a second method in the form of a Maximum Likelihood Estimator (MLE) to calculate the intrinsic scatter across our sample. This method produces larger uncertainties than our pairwise method, but it also has fewer sampling restrictions and can be applied to a larger set of OCs. The form of the MLE is shown in Equation <ref> below.

    lnL = ∏_i=1^∞1/√(2π)(σ _[X/Fe]^2 + e_i^2)^1/2exp(-(x_i-μ _[X/Fe])^2/2(σ _[X/Fe]^2 + e_i^2)).

In Equation <ref>, σ_ [X/Fe] is the intrinsic scatter being tested, and e_i is the uncertainty on the [X/Fe] measurement for the ith star in that cluster. We sample a narrow range of mean [X/Fe] (μ), where Δμ = 0.05 dex around the calculated mean of the cluster members, with an initial range of 0.1 dex for intrinsic scatter and spacing of 0.003 dex. We then do a second iteration with finer spacing in the intrinsic scatter dimension, with spacing on the order of 10^-4 dex, centered on the likeliest value from the coarser grid. An example of this is shown in Figure <ref>. We calculate a variance from the Fisher information matrix, and from that we calculate the uncertainty the intrinsic scatter. We apply the MLE method to all twenty-six clusters in our sample. Of these twenty-six clusters in our sample, ten
use the MLE intrinsic scatters for their final measurement.




  §.§.§ MLE Corrections for Small Samples


Based on tests with synthetic data, we find that the measured intrinsic scatter is unreliable when the number of cluster stars (N) is low, with a consistent systematic bias at N<9 in the measured scatter. This is due to the fact that the MLE is a biased estimator and at small sample sizes has a negative bias that causes it to underestimate the true parameter value, which can be accounted for with a multiplicative factor <cit.>: 


    Δσ_true = √(N/N-1)σ_MLE


This correction is applied to clusters with 6 ≤ N < 9 members. Below six members, the measured scatters are too unpredictably noisy. This small sample size correction results in an additional six clusters added to our sample, creating the final sample of 26 OCs.



  §.§.§ Systematics


As discussed in <cit.>, the methodology for measuring a star's [X/Fe] involves a multi-parameter fit that includes fitting the observed stellar spectrum to synthetic models. First a global fit to the spectrum is done to determine the best fit values for temperature, surface gravity, microturbulent velocity, and [M/H]. Holding these parameters constant, individual elemental abundances are extracted from narrow spectral windows. 

To test for any systematic trend between the global stellar parameters and elemental abundances, we quantify the slope of the uncalibrated log(g) vs. [X/Fe] within all the clusters in each pipeline, approximating it as a linear relationship, as seen in Figure <ref>. We exclude red clump stars as they are further along their evolutionary track than red giant stars, and potentially have slightly different surface abundances due to evolutionary effects or internal systematics. As a result, including them in the slope measurement could artificially drive any measurements of chemical homogeneity. However, these RC stars are adjusted afterwards along with the other giants in their cluster. This process ensures a uniform abundance correction across the entire cluster. From this, we find that the existing slopes in the cluster sample are nonzero, with a median slope of ∼0.02 dex/dex. 

This systematic bias is present in DR17, and in the Cannon and ASPCAP allStar files from IPL-3. However it is not present in the IPL-3 allStar release analyzed using the Payne. We adjust the measured [X/Fe] for each cluster star using the following equation:


    [X/Fe]_i, corrected = [X/Fe]_i - mlogg_i + ZP,


where each ith index is a star in a cluster and m is the best fit slope of [X/Fe] and log(g) within a specific cluster. We set the zero-point (ZP) of the cluster [X/Fe] after the correction using the abundances of the stars on the giant branch below the 
red clump to ensure that the median cluster [X/Fe] is reflective of its true value. The fitting uncertainties are propagated to uncertainties on the correction, which are then added in quadrature with the existing abundance uncertainties. 





§ RESULTS


Within one standard deviation, the only abundance that showed evidence of inhomogeneity, or consistent nonzero intrinsic scatter, across multiple clusters was [M/H]. This is because [M/H] has small uncertainties compared to other [X/Fe] measurements (∼0.008 dex as compared to ∼0.015 dex for the rest of the DR17 and IPL-3 abundances and ≥0.03 for the BAWLAS abundances). Given that the [M/H] uncertainties are smaller but not under-reported, as we verified in Section 2, it potentially implies that we could detect the existence of inhomogeneities that were then masked by larger uncertainties in the other abundances. 

However, within three standard deviations (a 99.7% confidence interval), none of the clusters show measurable inhomogeneities in any of the measured elements. Furthermore, the scale of the limits derived using the paired stars method are also on the scale of ∼0.001 dex in many elements. And using that method we find no measurement of scatter across in any of the elements and clusters. As a result, we are confident that in the majority of elements we can constrain the homogeneity of the OCs to less than 0.01 dex in a 99.7% confidence interval.

All the measured quantities for each [X/Fe] are presented in Table <ref>. We show that across DR17 and all MWM pipelines, we do not find any elements that show consistent chemical inhomogeneity across our cluster sample, and in Figure <ref> we show that we do not find any clusters with consistent scatter across their abundance samples. We only show the results from the IPL-3 ASPCAP data in the here because it includes values for weak lined elements. The literature comparison plots using DR17, IPL-3 Cannon, and IPL-3 Payne releases are comparable and are shown in Appendix <ref>. 




The upper limits on the intrinsic scatter measured in elements included in the BAWLAS catalog are higher than the limits on intrinsic scatter placed on the more well-measured elements, such as Mg or Ni. The reasons for this are twofold: Firstly, while many of the clusters studied did have enough stars to measure an intrinsic scatter, the number of stars that contained BAWLAS abundances within each cluster was smaller than the number of stars used to calculate α and iron-peak elements. Secondly, the associated uncertainties for the weak-lined elements were appreciably larger (0.03–0.08 dex) than the ones included in DR17 and MWM (0.01–0.04 dex). 





§ DISCUSSION




 §.§ Comparison to Milky Way Field Stars



To quantify the difference in chemical homogeneity between our OCs and the Milky Way field, we create a matched field star sample (MFS) that mirrors our existing cluster sample.

For each of the 15 clusters in our study with enough members to apply the pairwise method, we match each star within the cluster to a field star within two sigma[We tested the effects of using 1σ and 3σ to match stars and found no difference in our conclusions.] in Galactocentric radius, [M/H], [α/M], T_ eff, and logg. Here we consider the parameter uncertainties of both the cluster star as well as any candidate field stars. We use each star's [M/H] and [α/M] value from prior to the stellar parameter correction. We then measure the intrinsic scatter in each of the MFS samples, replicating the methodology outlined in Section <ref>. Finally, we compare the difference in intrinsic scatter between our MFS sample and the OC sample (Figure <ref>). 

We find that on average, across all abundances, the matched field star samples have +0.012_-0.01^+0.02 dex more intrinsic scatter than the open clusters in our sample. This is in strong agreement with <cit.>, which states that stars in the Milky Way are chemically similar (within ∼0.01–0.02 dex) when given a fixed Galactocentric radius, [M/H], and [α/Fe]. The median difference between OC intrinsic scatter and field star intrinsic scatter (Δσ_[X/Fe]) for each nucleosynthetic family is given below:

 -2pt
  
  * α-elements (Mg, Si, Ca, Ti, P, S): 0.002 dex
  
  * Iron-peak elements (Cr, Mn, Fe, Co, Ni, V): 0.012 dex
  
  * Odd-Z (Na, Al, K): 0.023 dex
  
  * Neutron-capture (Nd, Ce): 0.02 dex


Due to our selection criteria for the field star sample, we expect similar intrinsic scatter in the α, and iron-peak elements. 
Interestingly, two of the odd-Z elements, Al and K, both have nonzero scatter in the majority of our MFS samples despite being measured as homogeneous in our OCs. While in [Na/Fe] the distribution in Δσ_[X/Fe] is roughly symmetric, Na is an element with weak lines in APOGEE's wavelength range and was included in our sample with the BAWLAS catalog. As a result, not only does it have comparatively larger uncertainties than the other odd-Z elements, but it has poorer completion as well, as only a subset of high SNR stars in our study have BAWLAS abundances. This implies that odd-Z elements may be a useful tool in differentiating otherwise chemically similar populations. However, given that the distribution for all of these elements is consistent with zero in at least a subset of clusters and field star comparisons, more precise limits are needed to accurately test this 


Neutron capture elements also show slightly larger differences in scatter than their field star counterparts. <cit.> showed that neutron-capture elements have more discriminatory power in distinguishing “doppelganger stars”. What we find potentially corroborates that, but we also show that for co-eval and co-natal stars within an OC, the expectation of chemical homogeneity within neutron capture elements is broadly comparable to that of elements from other well-studied and well-measured nucleosynthetic families. In a field star sample with a high expectation of chemical similarity to an OC, as shown by the relative lack of difference in measured α-element intrinsic scatter, there is a 0.02 ± 0.02 dex difference in measured scatter for neutron capture elements. Due to this, it seems possible that neutron capture elements could be useful in distinguishing otherwise highly similar stellar populations but significantly more precise limits would be required to accomplish that goal.

In both cases of odd-Z and neutron capture elements, the differences between each OC and their respective MFS are consistent with zero within 3 uncertainties. More precise limits are required to make any conclusive statements on their effectiveness in distinguishing co-natal and co-eval stellar populations.



 §.§ Previous Studies of Chemical Homogeneity

There have been numerous studies focusing on measuring the chemical homogeneity of star clusters; however, most of these studies have been focused on larger and more complicated globular clusters. Within open clusters, different studies have found a wide range of limits on inhomogeneity. This is further complicated by the fact that not every study investigates the same set of abundances, nor is every analysis method comparable to one another. One of the differences between this study and many others is that their published limits on homogeneity are either 68% and 95% limits; the values we publish and show in Figures <ref> and Appendix <ref> are all 99.7% limits.

Given the number of studies done on M67, NGC 6791, and NGC 6819, they are shown in Figure <ref>, while the literature comparisons for the remaining clusters are located in Appendix <ref>. The figures comparing the intrinsic scatter measured using DR17, Cannon, and Payne abundances for all 26 OCs are also shown in Appendix <ref>.

We compare our results to four previous studies <cit.>. We find that within the α and iron peak elements, which are well-measured in APOGEE and MWM, the upper limits derived in this study are consistent with previous findings. In the weak-lined elements such as V, Cu, Ce, and Nd, there is more variance. But even within those elements, in many clusters we find comparable limits to previous works.

It is worth outlining the differences in the analyses and sample sets between these different studies. <cit.> is the most similar to ours in terms of sample size, analysis method, and dataset <cit.>. While they measure upper limits that are far less constraining than ours, they also measure lower limits that strongly imply the existence of true intrinsic scatter. However, the stellar parameter systematic that we found in DR17 and MWM is also present in DR16 but was unaccounted for in <cit.>'s final results. Therefore, it is possible that the measured scatters in <cit.> are impacted by a relationship between stellar parameters and abundance measurements. Our uncorrected σ_[X/Fe] values (not shown in this paper) are in strong agreement with the ones published in <cit.>, which lends evidence to this hypothesis.

<cit.> uses high resolution spectroscopy (R ∼ 50,000) to study two solar twins in M67. This, makes it less likely that their final abundances are driven by systematics due to stellar atmospheres or poor line fitting, as the stars are in very similar parts of parameter space. This method is similar, but not identical, to the pairwise method of deriving intrinsic scatter outlined in Section <ref>. They derive abundances for a total of 26 elements as well as [Fe/H], with an average measurement uncertainty of e_[X/Fe] ≤ 0.02 dex. As a result, within their sample they were able to very tightly constrain the difference in [X/Fe] between the two stars, showing that for elements with Z<30, there was an average difference of < 0.06 dex, and for the neutron capture elements there was an average difference of < 0.05 dex. It is encouraging that in the strong-lined elements included in our study, we place similar upper limits on abundance scatter.

<cit.> uses APOGEE DR13 <cit.> abundances and spectra to measure intrinsic dispersions within a set of seven open clusters, six of which are included in this study. However, their methodology in constraining scatter was noticeably different than this work. Using The Cannon <cit.>, they derive abundances and uncertainties in 20 different elements from APOGEE spectra and a training set of open cluster stars in APOGEE DR13. <cit.> notes that while their Cannon abundance measurements were broadly comparable to ASPCAP's, the uncertainties are between 20-50% smaller. Beyond that, using a chi-squared fit they determined that the uncertainties, a quadrature sum of formal abundance uncertainty and cross-validation uncertainty, are overestimated given the widths of calculated abundance distributions. Therefore, they derive a scaling factor to correct the uncertainties to match the theoretically predicated value; however, this methodology also introduces the risk of artificially down-scaling the measured limits on the intrinsic scatter. However, it is also worth noting that in the majority of elements, the values they publish are comparable to the ones derived in this work.

<cit.> studied the abundance spread of 15 different elements in three clusters, all of which are included in this study: M67, NGC 2420, and NGC 6819. The data came from APOGEE DR12 <cit.>, with cluster membership from <cit.>. However, the key difference between their study and this one is the methodology. <cit.> made the assumption that in the absence of any intrinsic chemical scatter, the main driver for variation in the photometric and spectroscopic attributes of OC stars is their mass, which can be modeled as a one dimensional sequence. They correct for any systematic variations in the spectra driven by mass, using T_ eff as a proxy. They then use detailed forward modeling of the spectra and Approximate Bayesian Computation to measure the intrinsic scatter as well as upper limits. Overall <cit.> find no indication of chemical inhomogeneity in any of the three clusters they studied; the upper limits they derived are largely in agreement with the ones calculated using the MLE method in this work. However, using the paired stars method we can constrain tighter upper limits in the majority of elements. Similar to studies discussed previously, the limits placed on the BAWLAS neutron capture elements by <cit.> are lower than ours.

Finally, <cit.> use spectroscopic data from APOGEE DR14 <cit.> to measure intrinsic scatter in M67, NGC 6791, and NGC 6819 in 15 different elemental abundances. The analysis method is very similar to the one outlined in <cit.>, though there are a few differences — notably, that they use DR14 instead of DR12, which includes several differences in the line lists <cit.>. Furthermore, unlike <cit.>, they use spectroscopic effective temperatures in their one-dimensional model as opposed to photometric effective temperatures. <cit.> found the clusters to be chemically homogeneous, placing upper limits comparable to this study across its sample of elements. While they measured fewer abundances than this work, in the three clusters studied the upper limits they derive are similar to <cit.>.

Thus, within the majority of the elements included in DR17 and IPL-3, such as the α-elements and iron-peak elements, the upper limits we calculate here are in agreement with what has been previously found across all the clusters studied. This lends credibility to the limits placed on the numerous clusters studied in this work that did not have previously derived limits in the literature. However, the limits derived in this study for the neutron capture elements are larger than what has been previously found in any of the literature. This is likely driven by the comparatively large uncertainties on the elements (0.03–0.08 dex).



§ CONCLUSION


The purpose of this study was to quantify the level of chemical homogeneity within the largest sample to date of Milky Way open clusters for a broad set of elements. Using SDSS-V Milky Way Mapper IPL-3 abundances and Gaia DR3 kinematics, we identify a sample of 26 open clusters with large enough membership to measure the intrinsic scatter in up to 20 elements. Using the abundance differences between paired stars along the HR diagram, as well as a Maximum Likelihood Estimator, we then measure the intrinsic scatter within each element for each cluster. We find the following:

 
  * We assemble a sample of 26 open clusters across a broad range of metallicity, age, mass, and galactic radii. Within a 99.7% confidence interval, we do not find any evidence of intrinsic scatter on the giant branch or in the red clump in any element across all the open clusters in our sample. 
 
 
  * Within the majority of abundances included in APOGEE DR17 and Milky Way Mapper IPL-3, we constrain the chemical homogeneity to ≤0.02 dex within a 99.7% confidence interval, and within ≤0.2 dex for the weak lined elements, such as those included in the DR17 BAWLAS catalog. Our limits are consistent with those in the literature for well-studied elements and clusters, and we add roughly a dozen clusters to this literature sample. Given the limited dataset in some of the elements, we recommend follow up measurements to better quantify their upper limits.
 
 
  * When compared to a sample of field stars with similar Galactocentric radii, [α/M], and [M/H], we find our OCs to be more chemically homogeneous, with an average difference of ∼0.012 dex between the two samples. This corroborates previous findings that the dimensionality of chemical enrichment of the Milky Way is low, and can likely be explained through a few processes. In the future this could be useful in placing constraints on radial mixing and azimuthal variations within the Milky Way.

 
  * We identify surface-gravity-dependent abundance shifts within APOGEE DR17 and Milky Way Mapper IPL-3 (corrected for in this analysis). This systematic needs to be accounted for in similar future work. We also find that the abundance uncertainties within both APOGEE and MWM are accurately estimated.

 
  * These findings have implications for attempts to implement chemical tagging, especially strong chemical tagging, specifically showing that within the light elements alone it is not possible to confidently separate field stars and co-natal stars given similar stellar parameters and Galactic radii. The tightest abundance variation constraints in OCs may also help set limits on the rate of binary interactions and planetary engulfment in different environments.

 





§ ACKNOWLEDGMENTS

We thank the anonymous referee for their helpful and insightful comments. This material is based upon work supported by the National Science Foundation under Grant No. AST-2206542. A.S. also acknowledges support from the University of Utah's J. Irvin and Norma K. Swigart First-Year Summer Graduate Research Fellowship. P.F. acknowledges support from NSF Grant AST-2206541, and K.C. acknowledges support from NSF Grant AST-2206543.

Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions. 
SDSS-IV acknowledges support and resources from the Center for High Performance Computing at the University of Utah. The SDSS website is <www.sdss4.org>.
SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, Center for Astrophysics | Harvard & Smithsonian, the Chilean Participation Group, the French Participation Group, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU) / University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatário Nacional / MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.

Funding for the Sloan Digital Sky Survey V has been provided by the Alfred P. Sloan Foundation, the Heising-Simons Foundation, the National Science Foundation, and the Participating Institutions. SDSS acknowledges support and resources from the Center for High-Performance Computing at the University of Utah. The SDSS web site is <www.sdss.org>. SDSS is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration, including the Carnegie Institution for Science, Chilean National Time Allocation Committee (CNTAC) ratified researchers, the Gotham Participation Group, Harvard University, Heidelberg University, The Johns Hopkins University, L’Ecole polytechnique fédérale de Lausanne (EPFL), Leibniz-Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Extraterrestrische Physik (MPE), Nanjing University, National Astronomical Observatories of China (NAOC), New Mexico State University, The Ohio State University, Pennsylvania State University, Smithsonian Astrophysical Observatory, Space Telescope Science Institute (STScI), the Stellar Astrophysics Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Illinois at Urbana-Champaign, University of Toronto, University of Utah, University of Virginia, and Yale University.

astropy <cit.>, scipy <cit.>,matplotlib <cit.>, numpy <cit.>, pandas <cit.>, and seaborn <cit.>.




§ 


Here we show the literature comparison plots for the well studied OCs in our sample, as discussed in Section <ref>.









§ OPEN CLUSTER SAMPLE MEMBERSHIP


Here we show the plots outlining the kinematic selection of our open cluster sample, as described in Section <ref>.



