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 interstellar medium pollution and gas mixing in progenitor molecular clouds. Using Sloan Digital Sky Survey SDSSV Milky Way Mapper and SDSSIV Apache Point Observatory Galaxy Evolution Experiment 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.
['Open star clusters', 'Stellar astronomy', 'Star clusters', 'Galactic archaeology', '1160', '1583', '1567', '2178', 'Astrophysics - Astrophysics of Galaxies', 'Astrophysics - Solar and Stellar Astrophysics']
A Comprehensive Study of Open Cluster Chemical Homogeneity Using APOGEE and Milky Way Mapper Abundances
{'A Comprehensive Study of Open Cluster Chemical Homogeneity using APOGEE and Milky Way Mapper Abundances': "Amaya Sinha, 1 Gail Zasowski, 1 Peter Frinchaboy, 2 Katia Cunha, 3, 4 Diogo Souto, 5 Jamie Tayar, 6 and Keivan Stassun 7  \n1 Department of Physics & Astronomy, University of Utah, 115 South 1400 East, Salt Lake City, UT 84112, USA 2 Department of Physics & Astronomy, Texas Christian University, Fort Worth, TX 76129, USA  \n3 Steward Observatory, University of Arizona, Tucson, AZ 85721, USA  \n4  \nObservat'orio Nacional, Rua General Jos'e Cristino, 77, 20921-400 S˜ao Crist'ov˜ao, Rio de Janeiro, RJ, Brazil 5 Departamento de F'ısica, Universidade Federal de Sergipe, Av. Marcelo Deda Chagas, Cep 49.107-230, S˜ao Crist'ov˜ao, SE, Brazil 6 Department of Astronomy, University of Florida, Bryant Space Science Center, Stadium Road, Gainesville, FL 32611, USA 7 Department of Physics and Astronomy, Vanderbilt University, VU Station 1807, Nashville, TN 37235, USA", '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.', '1. 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 (Matteucci & Francois 1989; Chiappini et al. 2001; Spitoni & Matteucci 2011; Weinberg et al. 2019). 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.  \nIn the age of large astronomical surveys such as GALAH (De Silva et al. 2015; Buder et al. 2021; Sheinis et al. 2015; Kos et al. 2017; Zwitter et al. 2021), LAMOST (Zhao et al. 2012), RAVE (Steinmetz et al. 2020a,b), APOGEE (Majewski et al. 2017), and Gaia  \n(Gaia Collaboration et al. 2016; Katz et al. 2023), 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 (Hawkins et al. 2020) to large populations of dispersed field stars (Ness et al. 2022).  \nStars 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 (Lada & Lada 2003). 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 (Freeman & Bland-Hawthorn 2002). This technique, known as chemical tagging, has been a strong motivator for studies of cluster chemistry.  \nWhile many studies support this assumption of OC chemical homogeneity (De Silva et al. 2007; Bovy 2016;  \nNess et al. 2018; Cheng et al. 2021), there has been work showing that at least some clusters are chemically inhomogeneous (Ness et al. 2018). For example, Geisler et al. (2012) 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 (Gratton et al. 2001). 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 Cunha et al. (2015), 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 5.  \nThere 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 (Smith & Lambert 1990). 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 (e.g.; Bressan et al. 2012; Lagarde et al. 2012; Choi et al. 2016a). NGC 6705 is an interesting OC regarding this effect, as it has been observed to also be enhanced in Na due to dredge up (Loaiza-Tacuri et al. 2023).  \nThe surface abundances of elements such as Mg can be affected by effects like mass transfer and atomic diffusion (Michaud et al. 2015; Souto et al. 2019). However, the latter only weakly impacts the upper giant branch. Lastly, mass transfer (e.g.; Milliman et al. 2015; Bastian et al. 2013; Abate et al. 2013) and pollution events such as planetary engulfment (e.g.; Pinsonneault et al. 2001; Laughlin & Adams 1997; Carlberg et al. 2010) can also alter a star's surface abundances.  \nHowever, 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 (Feng & Krumholz 2014); 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 Bhattarai et al. (2024) who find open clusters in FIRE-2 simulations (Hopkins et al. 2014, 2018; Wetzel et al. 2016) to have chemical scatter within 0.02 dex on average.  \nHowever, chemical inhomogeneity in real clouds could be due to effects not fully captured by simulations, re-  \nlated 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 (e.g.; Krumholz et al. 2019; Looney et al. 2006).  \nQuantifying 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.  \nThis 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 2 outlines the survey data, verification of the abundance uncertainties, and determination of the cluster membership. Section 3 details the methodology and calculation of the intrinsic scatter within each [X/Fe] across the final cluster sample. Section 4 presents the results of our work, and Section 5 compares our results to previous findings.", '2.1.1. 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 (SDSS-V; Kollmeier et al. 2017, J. Kollmeier, in prep ). 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 (Almeida et al. 2023).  \nSDSS-V/MWM uses two telescopes: the Sloan Foundation Telescope at APO (Gunn et al. 2006) and the duPont Telescope at LCO (Bowen & Vaughan 1973). Both are outfitted with nearly identical custom-built 300-fiber APOGEE spectrographs (Wilson et al. 2019), 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 (Pogge et al. 2020), which benefited from the adoption of a three-element corrector for the Sloan telescope at APO (Barkhouser et al. 2022).  \nWithin IPL-3, three different data pipelines were used to analyze the data taken from APO and LCO: The Payne (Ting et al. 2019), The Cannon (Ness et al. 2015), and the APOGEE Stellar Parameters and Abundances Pipeline (ASPCAP; Garc'ıa P'erez et al. 2016). Both The Payne and The Cannon are label-transfer methods that determine stellar labels from spectra after being  \ntrained on a set of spectra with known labels. They are differentiated by the fact that The Cannon is a datadriven 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.", '2.1.2. SDSS-IV/APOGEE': "In addition to IPL-3, we use abundances and RVs from the seventeenth and final data release (DR17; Abdurro'uf et al. 2022) of SDSS-IV's (Blanton et al. 2017) Apache Point Observatory Galaxy Evolution Experiment (APOGEE; Majewski et al. 2017), which contains over 700,000 stars. The initial targeting strategy for APOGEE-1 and APOGEE-2 are outlined in Zasowski et al. (2013) and Zasowski et al. (2017), respectively, and the final targeting for APOGEE-2 is outlined in Beaton et al. (2021) for APOGEE-2 north, and Santana et al. (2021) for APOGEE-2S.  \nThe details for the APOGEE data reduction pipeline are described in Nidever et al. (2015), and the details for ASPCAP are found in Garc'ıa P'erez et al. (2016). 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 Jonsson et al. (2020) and Gustafsson et al. (2008). The line lists used for DR17 are outlined in Smith et al. (2021), and the spectral fitting used for ASPCAP is described in Allende Prieto et al. (2006). The details describing the APOGEE spectral grids can be found in Osorio et al. (2020a), Osorio et al. (2020b), and Hubeny & Lanz (2017), and lastly, the details for Turbospectum can be found in Plez (2012a) and Plez (2012b).  \nWe also include abundances from the BACCHUS Analysis of Weak Lines in APOGEE Spectra (BAWLAS; Hayes et al. 2022), 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 STARFLAG or ASPCAPFLAG and analyzed using the BACCHUS code (Masseron et al. 2016), 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  \nthe BAWLAS VAC abundances and uncertainties for the following elements: Na, P, S, V, Cu, Ce, and Nd.  \nTwo separate uncertainties are reported for each BAWLASabundance measurement. One is the X FE ERR MEAS describing the measured uncertainty from the combined spectra using the same methodology as ASPCAP. The other is X FE ERR EMP , which describes the uncertainty derived from the spectral lines themselves. Here to remain consistent in our analysis we use X FE ERR MEAS as it is the closest to the uncertainty calculation method that we verify in Section 2.3.  \nThe calculation of the abundances in the BAWLAS catalog is outlined fully in Hayes et al. (2022). Each spectrum has an associated X SPECTRA FLAG , with values 0 , 1 , or 9 . Both 0 and 9 indicate either suspicion with the final fit or total failure, and only 1 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 X SPECTRA FLAG = 1.", '2.1.3. Quality Cuts': "We limit our sample in APOGEE DR17 to stars with VERR < 0.1 km s -1 and |VHELIO AVG| < 5000 km s -1 to ensure our stars have reliable radial velocities. We also restrict our sample to stars with VSCATTER < 1 km s -1 , to remove potential binaries within our sample (Badenes et al. 2018; Price-Whelan et al. 2020). Here VHELIO AVG refers to the average radial velocity derived from individual RVs that are drawn from cross-correlation of individual spectra with combined spectrum. VERR refers to the uncertainty on that radial velocity, and VSCATTER refers to the scatter of individual visit RVs around the average. To ensure the sample has reliable measurements we enforce a SNR > 50 . 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.  \nWe also exclude from our sample any stars that have [X/Fe] flags in more than 2 elements. Lastly, we only sample from stars with LOGG ≤ 3.5 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 VSCATTER which is not included in IPL-3.  \nIn DR17's allStar file, we enforce quality cuts on the STARFLAG and ASPCAPFLAG 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  \nASPCAPFLAG column we enforced the following requirements:  \n- 1. BITMASK 23; STAR BAD == 0; BAD overall for star: set if any of TEFF , LOGG , CHI2 , COLORTE , ROTATION , SN error are set, or any parameter is near grid edge ( GRIDEDGE BAD is set in any PARAMFLAG )\n- 2. BITMASK 19; METALS BAD ==0; FERRE failed to return a value for metals.\n- 3. BITMASK 20; ALPHAFE BAD == 0; Elemental abundance from window differs > 0.5 dex from parameter abundance for [ α /Fe].  \nand within the STARFLAG we made the following cuts:  \n- 1. BITMASK 2; BRIGHT NEIGHBOR == 0; Star has neighbor more than 10 times brighter.\n- 2. BITMASK 3; VERY BRIGHT NEIGHBOR == 0; Star has neighbor more than 100 times brighter.\n- 3. BITMASK 17; SUSPECT BROAD LINES ==0; WARNING: cross-correlation peak with template significantly broader than autocorrelation of template.  \nThe IPL-3 ASPCAP allStar file does not have a published ASPCAPFLAG column. As such, we require that all the stars in our sample have flag bad == False and flag warn == False . For all three IPL-3 pipelines we use we enforce a result flags == 0 requirement. Where possible we enforce a x h flag == 0 requirement within each pipeline. Within the three IPL-3 allStar files we use, we enforce these quality cuts on the spectrum flag column, which has bits that correspond to DR17's STARFLAG column.", '2.2. Cluster Membership': "We start from the catalog of cluster members published in Cantat-Gaudin et al. (2018), hereafter CG18, which contains membership information for over 200,000 stars across ∼ 2000 OCs using Gaia DR2 (Gaia Collaboration et al. 2018). Due to the more recent availability of Gaia DR3 (Gaia Collaboration et al. 2021, 2023), we first match the stars identified as cluster members in CG18 to their DR3 kinematics and positions. CG18 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 CG18 within three cluster radii  \nof their cluster centers. Using each cluster's distribution in radial velocity, we find that a minimum probability cut of P ≥ 0.5 in the CG18 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 CG18, only 145 clusters have any members within both DR17 and IPL3.", '2.2.1. 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 1, 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 B.  \nWe 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 (Dotter 2016; Choi et al. 2016b), generated using ages from Cantat-Gaudin et al. (2020), 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 3.2.1.", '2.2.2. Final Cluster Sample': "The distribution of the final cluster sample in radius, age, and [Mg/Fe] versus [Fe/H] can be seen in Figure 2. We calculate the positional, kinematic, and orbital information for each cluster in our final cluster sample. Using Astropy SkyCoords (Astropy Collaboration et al. 2013, 2018, 2022) we calculate the Cartesian X, Y, Z galactocentric coordinates for all the clusters in our sample, as well as the galactocentric radius.  \nWe integrate each cluster's orbit to measure its Z max , eccentricity, and guiding radius using Galpy, with the MWPotential2014 gravitational potential (Bovy 2015; Mackereth & Bovy 2018) 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 Carrera et al. (2022), as we  \nFigure 1. Cluster membership plot for M67. Top Left: The distribution of stars on the sky in RA and DEC, with cluster members colored in blue, cluster member candidates that did not pass Section 2.2.1's kinematic and parallax cuts in red, and field stars in grey. Top Right: The KDE distribution of radial velocities of the cluster members compared to nearby field stars. Middle Left: The distribution of the combined proper motion of the cluster members compared to nearby field stars. Middle Right: The distribution of parallaxes of the cluster members compared to nearby field stars. Bottom Left: The distribution of [Fe/H] of the cluster members compared to nearby field stars; note that this was not used to prune membership. Bottom Right: The Kiel diagram of cluster members compared to surrounding stars. We also plot the MIST isochrone (Dotter 2016; Choi et al. 2016b) for this cluster's age and metallicity in black.  \n<!-- image -->  \ndo not use these parameters in the final science results, we do not include spiral arms in our potential.  \nWe 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 Weinberg et al. (2019), we quantify the ratio of nucleosynthetic enrichment within each cluster from CCSNe and Type Ia supernovae.  \nLastly, we include age estimates on all our clusters from Cantat-Gaudin et al. (2020). These ages were derived using an artificial neural network trained on a sample of reference clusters. For further description see  \nCantat-Gaudin et al. (2020). At this stage we also flag red clump stars by eye within each cluster. All the cluster parameters are included in Table 1, and will be included as a machine readable table.", '2.3. Verification of Abundances & Uncertainties': "The values of intrinsic scatter in OC abundances are on the order of ∼ 0.01 dex (e.g., Bovy 2016; Liu et al. 2016; Ness et al. 2018; Poovelil et al. 2020), 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 Poovelil et al. (2020). 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 1 from Poovelil et al. (2020).  \ne [ X/Fe ] ,k = √ π 2 median( | [X / Fe] j -[X / Fe] i ] | ) , (1)  \nwhere [ X/Fe ] j and [ X/Fe ] i are the abundance measurements from the same star's i th and j th visits. e [ X/Fe ] ,k is the resulting [X/Fe] uncertainty after median stacking the pairwise measurements in the k th bin.  \nWe group stars by SNR in bins spanning 50-70, 70-100, 100-130, 130-200, and greater than 200, as shown in Fig 3. 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.  \nWithin 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.", '3. MEASURING CHEMICAL SCATTER': "3.1. Paired Stars Method  \n<!-- image -->  \nFigure 2. Left : A plot of the [Mg/Fe] versus [Fe/H] plane  \n<!-- image -->  \nwhere the full giant star sample in APOGEE DR17 (subject to the same quality cuts as in Section 2.1.3) is shown in black, and the clusters are colored by log(age). Right : The guiding radius and galactocentric height of all the clusters in our sample, colored by number of member stars as described in Section 2.2.1.Figure 3. Top : [Mg/Fe] uncertainties calculated using stars with multiple visits. Middle : [Mg/Fe] uncertainties calculated by ASPCAP for DR17. Bottom : The subtracted difference between the uncertainties calculated using multiple visits and the uncertainties from ASPCAP.  \n<!-- image -->  \nTable 1. List of columns in the table of cluster parameters.  \nOur 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 3.2.2 is on the order ∼ 0.001 dex.  \nWe measure the intrinsic dispersion within each pair using Equation 2, 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.  \nσ = √ π 2 | ∆[ X/Fe ] | 2 -( e 2 1 + e 2 2 )) 2 (2)  \nWe 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 subsample, 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.  \nTo 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.  \nWe 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.  \nDue to the restrictions outlined above regarding the separation of the pairs in T eff and log g , as well as the minimum number of pairs required, this method can only be applied to 15 of the 26 clusters studied in this  \npaper. However due to the significantly higher precision of these values as compared to those derived using the method outlined in Section 3.2, within these 15 clusters we only publish results from this method.", '3.2. 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 3 below.  \nlnL = ∞ ∏ i=1 1 √ 2 π ( σ 2 [X / Fe] +e 2 i ) 1 / 2 exp( -(x i -µ [X / Fe] ) 2 2( σ 2 [X / Fe] +e 2 i ) ) . (3)  \nIn Equation 3, σ [X / Fe] is the intrinsic scatter being tested, and e i is the uncertainty on the [X/Fe] measurement for the i th 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 4. 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.', '3.2.1. 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 (Thorley & Schwarz 2018):  \n∆ σ true = √ N N -1 σ MLE (4)  \nThis 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.  \nFigure 4. Outside Red: The initial likelihood surface over which we initially test to determine a best guess for σ [ X/Fe ] . Within Red: Using the findings from the initial run, we more finely sample the likelihood surface in order to measure the final intrinsic scatter and uncertainty.  \n<!-- image -->', '3.2.2. Systematics': "As discussed in Garc'ıa P'erez et al. (2016), the methodology for measuring a star's [X/Fe] involves a multiparameter 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.  \nTo 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 5. 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.  \nThis systematic bias is present in DR17, and in the Cannon and ASPCAP allStar files from IPL-3. However  \nit 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:  \n[ X/Fe ] i, corrected = [ X/Fe ] i -m log g i +ZP , (5)  \nwhere each i th 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.  \nFigure 5. Uncalibrated surface gravity versus [Mg/Fe] for M67. The [Mg/Fe] values from DR17 are shown in black, and have a clear linear relationship with log(g). The adjusted abundances used in the calculation of intrinsic scatter are shown in red.  \n<!-- image -->", '4. 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  \npotentially implies that we could detect the existence of inhomogeneities that were then masked by larger uncertainties in the other abundances.  \nHowever, 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.  \nAll the measured quantities for each [X/Fe] are presented in Table 2. 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 6 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 A.  \nThe 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.010.04 dex).', '5.1. 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.  \nFor 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 1  \n1 We tested the effects of using 1 σ and 3 σ to match stars and found no difference in our conclusions.  \nin Galactocentric radius, [M/H], [ α /M], T eff , and log g . 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 3.1. Finally, we compare the difference in intrinsic scatter between our MFS sample and the OC sample (Figure 7).  \nWe find that on average, across all abundances, the matched field star samples have +0.012 +0 . 02 -0 . 01 dex more intrinsic scatter than the open clusters in our sample. This is in strong agreement with Ness et al. (2022), 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:  \n- · α -elements (Mg, Si, Ca, Ti, P, S): 0.002 dex\n- · Iron-peak elements (Cr, Mn, Fe, Co, Ni, V): 0.012 dex\n- · Odd-Z (Na, Al, K): 0.023 dex\n- · Neutron-capture (Nd, Ce): 0.02 dex  \nDue to our selection criteria for the field star sample, we expect similar intrinsic scatter in the α , and ironpeak 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  \nNeutron capture elements also show slightly larger differences in scatter than their field star counterparts. Manea et al. (2023) 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 conatal stars within an OC, the expectation of chemical homogeneity within neutron capture elements is broadly  \nTable 2. List of [X/Fe] columns in the table of cluster parameters for IPL-3 ASPCAP. All columns have units of dex, and all X/Fe measurements are scaled to solar values. The abundances measured in these columns are Mg, Ca, Si, S, P, Ni, Cr, Ti, V, Mn, Co, Na, K, Al, V, Mn, Co, Cu, CNO, Nd, and Ce. In addition we provide columns for [Fe/H] and [M/H]. This table will be made available in machine readable format.  \ncomparable 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.  \nIn 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.", '5.2. 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 6 and Appendix A are all 99.7% limits.  \nGiven the number of studies done on M67, NGC 6791, and NGC 6819, they are shown in Figure 6, while the literature comparisons for the remaining clusters are located in Appendix A. The figures comparing the intrinsic scatter measured using DR17, Cannon , and Payne  \nabundances for all 26 OCs are also shown in Appendix A.  \nWe compare our results to four previous studies (Bovy 2016; Liu et al. 2016; Ness et al. 2018; Poovelil et al. 2020). 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.  \nIt is worth outlining the differences in the analyses and sample sets between these different studies. Poovelil et al. (2020) is the most similar to ours in terms of sample size, analysis method, and dataset (using APOGEE DR16; Ahumada et al. 2020). 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 Poovelil et al. (2020)'s final results. Therefore, it is possible that the measured scatters in Poovelil et al. (2020) 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 Poovelil et al. (2020), which lends evidence to this hypothesis.  \nLiu et al. (2016) 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 3.1. 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  \n<!-- image -->  \n8  \nFigure 6. The 99.7% upper limits using abundances from IPL-3 ASPCAP in M67, NGC 6791, and NGC 6819. The median 99.7 percentile upper limit for each cluster is shown as a magenta horizontal line. We compare our findings to upper limits and measurements from previous literature. The upper limits are broadly in agreement with other studies for well-measured elements in APOGEE and Milky Way Mapper.Figure 7. A violin plot comparing the difference in element intrinsic scatter between open clusters and matched samples of field stars as shown in Section 5.1. Here the elements that belong to the α , CNO, and iron-peak families are shaded in grey as we expect them to be similar due to our selection criteria. Within these elements P and V have larger limits due to them being weak lined elements from the BAWLAS catalog. While broadly we find that the field star samples have more scatter than their OC counterpart, we don't find any conclusive difference in any element between the two samples.  \n<!-- image -->  \ntheir 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.  \nNess et al. (2018) uses APOGEE DR13 (Albareti et al. 2017) 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 (Ness et al. 2015), they derive abundances and uncertainties in 20 different elements from APOGEE spectra and a training set of open cluster stars in APOGEE DR13. Ness et al. (2018) 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.  \nBovy (2016) 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 (Alam et al. 2015), with cluster membership from M'esz'aros et al. (2013). However, the key difference between their study and this one is the methodology. Bovy (2016) 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 Bovy (2016) 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  \nlimits in the majority of elements. Similar to studies discussed previously, the limits placed on the BAWLAS neutron capture elements by Bovy (2016) are lower than ours.  \nFinally, Cheng et al. (2021) use spectroscopic data from APOGEE DR14 (Abolfathi et al. 2018) 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 Bovy (2016), though there are a few differences - notably, that they use DR14 instead of DR12, which includes several differences in the line lists (detailed in Holtzman et al. 2018). Furthermore, unlike Bovy (2016), they use spectroscopic effective temperatures in their one-dimensional model as opposed to photometric effective temperatures. Cheng et al. (2021) 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 Bovy (2016).  \nThus, within the majority of the elements included in DR17 and IPL-3, such as the α -elements and ironpeak 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).", '6. 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:  \n- 1. 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.  \n- 2. Within the majority of abundances included in APOGEE DR17 and Milky Way Mapper IPL3, 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.\n- 3. 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.\n- 4. 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.\n- 5. 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.', '7. 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 AST2206541, and K.C. acknowledges support from NSF Grant AST-2206543.  \nFunding 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 fur Astrophysik Potsdam (AIP), Max-Planck-Institut fur Astronomie (MPIA Heidelberg), Max-Planck-Institut fur Astrophysik (MPA Garching), Max-Planck-Institut fur Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observat'ario Nacional / MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Aut'onoma de M'exico, 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.  \nFunding 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'ed'erale de Lausanne (EPFL), Leibniz-Institut fur Astrophysik Potsdam (AIP), Max-Planck-Institut fur Astronomie (MPIA Heidelberg), Max-Planck-Institut fur Extraterrestrische Physik (MPE), Nanjing University, National Astronomical Observatories of China (NAOC), New Mexico State University, The Ohio State University, Pennsylvania  \nState University, Smithsonian Astrophysical Observatory, Space Telescope Science Institute (STScI), the Stellar Astrophysics Participation Group, Universidad Nacional Aut'onoma de M'exico, University of Arizona, University of Colorado Boulder, University of Illinois at  \nUrbana-Champaign, University of Toronto, University of Utah, University of Virginia, and Yale University.  \nSoftware: astropy (Astropy Collaboration et al. 2013, 2018,2022),scipy(Virtanenetal.2020),matplotlib(Hunter 2007), numpy (Harris et al. 2020), pandas (pandas development team 2020; Wes McKinney 2010), and seaborn (Waskom2021).", 'REFERENCES': "- -. 2012b, Turbospectrum: Code for spectral synthesis, Astrophysics Source Code Library, record ascl:1205.004. http://ascl.net/1205.004  \nPogge, R. W., Derwent, M. A., O'Brien, T. P., et al. 2020, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 11447, Ground-based and Airborne Instrumentation for Astronomy VIII, ed. C. J. Evans, J. J. Bryant, & K. Motohara, 1144781,  \ndoi: 10.1117/12.2561113  \n- Poovelil, V. J., Zasowski, G., Hasselquist, S., et al. 2020, ApJ, 903, 55, doi: 10.3847/1538-4357/abb93e\n- Price-Whelan, A. M., Hogg, D. W., Rix, H.-W., et al. 2020, ApJ, 895, 2, doi: 10.3847/1538-4357/ab8acc\n- Santana, F. A., Beaton, R. L., Covey, K. R., et al. 2021, AJ, 162, 303, doi: 10.3847/1538-3881/ac2cbc\n- Sheinis, A., Anguiano, B., Asplund, M., et al. 2015, Journal of Astronomical Telescopes, Instruments, and Systems, 1, 035002, doi: 10.1117/1.JATIS.1.3.035002\n- Smith, V. V., & Lambert, D. L. 1990, ApJS, 72, 387, doi: 10.1086/191421\n- Smith, V. V., Bizyaev, D., Cunha, K., et al. 2021, AJ, 161, 254, doi: 10.3847/1538-3881/abefdc\n- Souto, D., Allende Prieto, C., Cunha, K., et al. 2019, ApJ, 874, 97, doi: 10.3847/1538-4357/ab0b43  \nSpitoni, E., & Matteucci, F. 2011, A&A, 531, A72,  \ndoi: 10.1051/0004-6361/201015749  \n- Steinmetz, M., Matijeviˇc, G., Enke, H., et al. 2020a, AJ, 160, 82, doi: 10.3847/1538-3881/ab9ab9  \n- Steinmetz, M., Guiglion, G., McMillan, P. J., et al. 2020b,\n- AJ, 160, 83, doi: 10.3847/1538-3881/ab9ab8 Thorley, J., & Schwarz, C. 2018, Journal of Open Source Software, 3, 1082, doi: 10.21105/joss.01082 Ting, Y.-S., Conroy, C., Rix, H.-W., & Cargile, P. 2019, ApJ, 879, 69, doi: 10.3847/1538-4357/ab2331 Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2 Waskom, M. L. 2021, Journal of Open Source Software, 6, 3021, doi: 10.21105/joss.03021 Weinberg, D. H., Holtzman, J. A., Hasselquist, S., et al. 2019, ApJ, 874, 102, doi: 10.3847/1538-4357/ab07c7 Wes McKinney. 2010, in Proceedings of the 9th Python in Science Conference, ed. St'efan van der Walt & Jarrod Millman, 56 - 61, doi: 10.25080/Majora-92bf1922-00a Wetzel, A. R., Hopkins, P. F., Kim, J.-h., et al. 2016, ApJL, 827, L23, doi: 10.3847/2041-8205/827/2/L23 Wilson, J. C., Hearty, F. R., Skrutskie, M. F., et al. 2019, PASP, 131, 055001, doi: 10.1088/1538-3873/ab0075 Zasowski, G., Johnson, J. A., Frinchaboy, P. M., et al. 2013, AJ, 146, 81, doi: 10.1088/0004-6256/146/4/81 Zasowski, G., Cohen, R. E., Chojnowski, S. D., et al. 2017, AJ, 154, 198, doi: 10.3847/1538-3881/aa8df9 Zhao, G., Zhao, Y., Chu, Y., Jing, Y., & Deng, L. 2012, arXiv e-prints, arXiv:1206.3569, doi: 10.48550/arXiv.1206.3569  \nZwitter, T., Kos, J., Buder, S., et al. 2021, MNRAS, 508, 4202, doi: 10.1093/mnras/stab2673", 'APPENDIX': "A.  \nHere we show the literature comparison plots for the well studied OCs in our sample, as discussed in Section 5.  \nFigure 8. The literature comparison plots for the well studied OCs in our sample using IPL-3 ASPCAP (colored in magenta). The horizontal line indicates the median 99.7% confidence limit for that cluster's homogeneity.  \n<!-- image -->  \nFigure 9. The literature comparison plots for the well studied OCs in our sample using IPL-3 Cannon (colored in blue). The horizontal line indicates the median 99.7% confidence limit for that cluster's homogeneity.  \n<!-- image -->  \nFigure 10. The literature comparison plots for the well studied OCs in our sample using IPL-3 Payne (colored in red). The horizontal line indicates the median 99.7% confidence limit for that cluster's homogeneity.  \n<!-- image -->  \nFigure 11. The literature comparison plots for the well studied OCs in our sample using APOGEE DR17 (colored in black). The horizontal line indicates the median 99.7% confidence limit for that cluster's homogeneity.  \n<!-- image -->", 'B. OPEN CLUSTER SAMPLE MEMBERSHIP': 'Here we show the plots outlining the kinematic selection of our open cluster sample, as described in Section 2.2.1.  \nFigure 12. These figures follow the same layout as Figure 1 in Section 2.2.1.  \n<!-- image -->  \nFigure 13. These figures follow the same layout as Figure 1 in Section 2.2.1.  \n<!-- image -->  \nFigure 14. These figures follow the same layout as Figure 1 in Section 2.2.1.  \n<!-- image -->'}