Distant Relatives: The Chemical Homogeneity of Comoving Pairs Identified in Gaia

Stars that are born from the same gas, i.e., conatal, are critical in astronomy. The conatal nature of open clusters and wide binaries have made them indispensable laboratories for testing and improving our understanding of various areas of Galactic and stellar astrophysics. Both open clusters and wide binaries have been used to evaluate the current feasibility of chemical tagging (De Silva et al. 2007; Ting et al. 2012; Ness et al. 2018; Andrews et al. 2019; Hawkins et al. 2020b). Open clusters are used to calibrate stellar parameter and chemical abundance pipelines for large surveys (e.g., García Pérez et al. 2016). Wide binaries can be used to calibrate the metallicity of M-dwarfs (e.g., Lépine & Bongiorno 2007; Rojas-Ayala et al. 2010; Montes et al. 2018), constrain the age-magnetic activity relation (e.g., Garcés et al. 2011; Booth et al. 2017), the age–metallicity relation (e.g., Rebassa-Mansergas et al. 2016), and the initial-final mass relation (e.g., Zhao et al. 2012; Andrews et al. 2015). Besides, we can also use wide binaries to study exoplanet engulfment (e.g., Meléndez et al. 2017; Oh et al.2018).

Gravitationally bound binaries favor a conatal origin. For separations within a few thousand astronomical units, binaries are mostly formed via turbulent core fragmentation (Offner et al. 2010; Lee et al. 2017). At separations up to ∼10⁵ au, wide binaries can be formed from dynamical unfolding of triple systems (Reipurth & Mikkola 2012), the evaporation of star clusters (Kouwenhoven et al. 2010; Moeckel & Clarke 2011), and chance gravitational capture of prestellar cores (Tokovinin 2017). The predicted conatal nature of wide binaries motivated studies of their chemistry (Andrews et al. 2019; Ramírez et al. 2019; Hawkins et al. 2020b). They have found wide binaries to be chemically homogeneous to ∼0.02–0.1 dex in Δ[Fe/H] ⁸ Nonetheless, the studies of chemical homogeneity for binaries remain rather limited (${ \mathcal O }(10)$ pairs) as large surveys often do not observe both stars due to the subsampling, and we have to resort to targeted studies.

Comoving pairs, on the other hand, are even less studied than wide binaries. The investigation of detailed chemistry for comoving pairs remains largely nonexistent for separations beyond a few parsecs. Comoving pairs separated by ≳10⁵ au are no longer bound (Jiang & Tremaine 2010). As outlined in Oh et al. (2017), these unbound systems could arise from disrupted/dissolved conatal systems (open clusters, stellar binaries) or unrelated systems (chance alignments, gravitational resonance). Recently, simulations from Kamdar et al. (2019a) argued comoving pairs have a high probability of being conatal provided that their 3D velocity difference (Δv_3D) is below 2 km s⁻¹ at separations up to ${ \mathcal O }({10}^{6})$ au. Their simulations predicted that pairs with separations ∼10⁶ au are ∼80% conatal. If true, this could dramatically expand the number of systems that could be conatal, which would improve stellar atmospheric models, studies of exoplanet engulfment, and calibrate surveys, among other applications. With the precise astrometry provided by the recent data release of Gaia (Gaia Collaboration et al. 2018, 2021) sufficient comoving systems can readily be located, and this study is set up to test this proposition—are comoving stars with large separations conatal?

A key characteristic of conatal stars is their chemical compositions homogeneity; conatal stars are expected to have similar initial chemical compositions as they are usually formed from well-mixed interstellar medium (ISM) gas (Feng & Krumholz 2014). Nonetheless, measured surface abundances can be changed by internal processes (atomic diffusion, rotational mixing, dredge-up, gravitational settling, radiative levitation; Schuler et al. 2011; Dotter et al. 2017), model limitations (incomplete laboratory data, NLTE, 3D (Ruchti et al. 2013; Heiter et al. 2015; Nissen & Gustafsson 2018; Jofré et al. 2019), and methodological limitations (Jofré et al. 2014, 2019). All these effects could complicate any interpretation of the chemical homogeneity of the two stars.

To distinguish model systematic from astrophysical processes, in this study, we will appeal to differential analysis of stellar twins, i.e., stellar pairs with similar stellar parameters. Differential abundances can remove many of these effects provided the stars are close in stellar parameters because the systematic effects in the star of interest and reference star largely cancel out (see Gray 2005; Nissen & Gustafsson 2018; Jofré et al. 2019).

In this study, we measure abundance differences in 24 elements for 33 comoving pairs with Δv_3D < 2 km s⁻¹ and separations between ∼10² and 10⁷ au. With access to homogeneous high-resolution, high signal-to-noise data, we uniquely positioned to explore the detailed abundances from the wide binaries and the comoving stars and the chemical homogeneity of the two components. In Section 2, we describe the observations. Section 3 details the methods for estimating stellar parameters and abundances. Our results are presented in Section 4, in which we explore the conatal fraction of the comoving pairs. These results are discussed in Section 5 and summarized in Section 6.

Our target selection focused on comoving pairs of F and G dwarfs within 300 pc of the solar neighborhood. F and G dwarfs were selected because they are both luminous, and their spectral models are best characterized. We create the initial list of sources with the following ADQL query:

• SELECT *
• FROM gaiadr2.gaia_source
• WHERE ra > 165
• OR ra < 15
• AND dec < 25
• AND bp_rp BETWEEN 0.5 AND 1.5
• AND 1000/parallax < 300
• AND parallax_over_error > 5
• AND phot_g_mean_mag IS NOT NULL
• AND ((phot_g_mean_mag---5*LOG10(1000/parallax) + 5) BETWEEN 2 AND 7)
• AND radial_velocity IS NOT NULL

Our range of bp_rp was selected to exclude hot stars as some of them could be rapid rotators with broad spectral features that complicate spectral analysis. Next, to ensure that our systems are genuine (isolated) comoving stars, we filter out members of open clusters and stellar associations with a friend-of-friend (FoF) connectivity cut following Kamdar et al. 2021, removing any pairs that belong to an aggregate with more than four FoF companions (with separations <10 pc, Δv_3D < 5 km s⁻¹). We verify the efficacy of this filtering by crossmatching our sample with the open clusters and comoving groups catalogs from Kounkel & Covey (2019) and Cantat-Gaudin & Anders (2020), finding no overlap. From this list of possible sources, we calculate the 3D velocity separation Δv_3D and spatial separation using galpy (Bovy 2015) using astrometry and radial velocity (RV) measurements from Gaia DR2. The target selection prioritized stars in similar locations on the color–magnitude diagram (CMD) to reduce the influence of systematic errors on the derived elemental abundances. Together, we consider the comoving pairs that satisfy the following criteria:

• 1.  
  ∣ΔM_G∣ < 0.5 mag
• 2.  
  ∣Δ(G_BP − G_RP)∣ < 0.1 mag
• 3.  
  Δv_3D < 2 km s⁻¹
• 4.  
  3D separation < 300 pc

The remaining comoving pairs were binned in log separation and sorted by apparent magnitude. The final sample was selected by taking the most luminous ∼10 pairs from each log separation bin.

Throughout this study, we will use the term comoving stars to signify both bound and unbound pairs of stellar companions. In particular, for consistency, we will use the words close comoving pairs and wide binaries interchangeably to refer to (mostly bound) wide binaries and far comoving pairs to be the primarily unbound systems. For simplicity, we define close comoving pairs as those with 3D spatial separation of 1 pc (i.e., 2 × 10⁵ au), and the others far comoving pairs. This choice is partly motivated by the fact that in simulations from Kamdar et al. (2019b), the minimum separation of unbound comoving pairs have a separation of 1 pc. This is also consistent with the previous study from Jiang & Tremaine (2010).

We observed 33 pairs of comoving stars using the MIKE spectrograph on board the Magellan/Clay telescope (Bernstein et al. 2003) from 2019 June 13–16. These 33 pairs constitute our main sample. We also observe two additional pairs with Δv_3D > 2 km s⁻¹ as a control sample. But unless otherwise stated, we will only refer to the main sample throughout this study. Additionally, we also observed four Gaia benchmark stars (Jofré et al. 2014; Blanco-Cuaresma et al. 2014b; Heiter et al. 2015) dispersed over the sampled CMD to improve the data reduction. The benchmark stars used were 18 Sco, bet Vir, HD 140283, and mu Ara. The instrument employs a blue and red spectrograph to cover 3350–5000 Å and 4900–9500 Å respectively. We used the 0.5' slit with 2 × 1 binning, which gave the blue and red spectrographs' typical resolving powers of 50,000 and 40,000, respectively. The median signal-to-noise ratio (S/N) per pixel for the blue and red chips were 121 and 185. The observational details are given in Table 1.

Table 2. The 3D Separations (S) and Velocity Differences (Δv_3D) between Comoving Pairs

Gaia eDR3 ID for	Gaia eDR3 ID for	Identifier	SeDR3	${\sigma }_{{{\rm{S}}}_{\mathrm{eDR}3}}$	Δv3D	${\sigma }_{{\rm{\Delta }}{v}_{3{\rm{D}}}}$	Sgeometric	${\sigma }_{{\mathrm{lower}}_{{{\rm{S}}}_{\mathrm{geometric}}}}$	${\sigma }_{{\mathrm{upper}}_{{{\rm{S}}}_{\mathrm{geometric}}}}$
Component A	Component B		(au)	(au)	(km s−1)	(km s−1)	(au)	(au)	(au)
6160638833832941312	3471119180522314624	CM00	1.1 × 106	2.7 × 104	0.86	0.05
3613454637229633664	3613454637229634176	CM01	9.0 × 104	1.8 × 105	1.02	0.08	8.5 × 103	8.3 × 102	3.7 × 103
3639520621950395904	3639520621950395776	CM02	2.2 × 104	9.3 × 103	1.29	0.04	2.8 × 102	27	1.2 × 102
5897704951769034368	5897785667103983616	CM03	1.5 × 106	1.0 × 105	0.93	0.05
6224633983987510528	6224633983987511552	CM04	1.1 × 104	1.6 × 104	1.71	0.08	4.2 × 102	41	1.9 × 102
6004256909931771136	6004256802548180736	CM05	1.2 × 103	8.6 × 103	0.82	0.06	8.2 × 102	80	3.6 × 102
5877155289930033024	5877155289930029184	CM06	6.5 × 104	6.7 × 104	0.70	0.06	3.7 × 103	3.6 × 102	1.6 × 103
5798991008295109120	5798991008295120896	CM07	3.2 × 105	7.1 × 104	1.43	0.07
6214634887804201728	6214634887804233088	CM08	6.4 × 103	4.8 × 104	0.58	0.04	5.0 × 102	49	2.2 × 102
6208919660724640896	6208919660723526272	CM09	8.1 × 104	3.2 × 104	0.79	0.06	8.3 × 102	81	3.6 × 102

Note. A full machine-readable version is available online. The subscript eDR3 indicates separations calculated solely from Gaia eDR3 astrometry and the spectroscopic radial velocity (RV). As discussed in Section 2.1, the separations calculated this way may overestimate the true separation between the close comoving components. For the subset of comoving pairs suspected to be wide binaries, we invoke a geometric prior based to infer the 3D separations based on the better-measured 2D projected separations (Appendix A). We only perform this correction for the close comoving pairs. These separations are denoted with the subscript geometric. The uncertainty for the geometric separations is asymmetric, so the upper and lower error bars (16 and 84 percentiles) for this separation are listed separately. The results in this paper assumes S_geometric for close comoving pairs (i.e., S_eDR3 < 2 × 10⁵ au) and S_eDR3 for the far comoving pairs.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

Our sample size of 33 pairs is comparable to previous detailed chemical studies of wide binaries (Andrews et al. 2019; Ramírez et al. 2019; Hawkins et al. 2020b), but expands the range of separations by two orders of magnitude. Previous studies mostly restrict their sample to ≲10⁵ au, whereas in this study we have comoving pairs with separations up to 10⁷ au. In the top panel of Figure 1 we show the CMD of our comoving targets (white circles) with a random sample of 200,000 stars from Gaia DR2 plotted in the background for reference.

Zoom In Zoom Out Reset image size

Among our sample, we found that, for two particular pairs, each of which had a component that exhibit large deviations in the RV measurements from Gaia eDR3 ⁹ and our RV derived with the MIKE spectra. Both systems were flagged as outliers based on the interquartile range (IQR) test. ¹⁰ This deviation could indicate an unresolved companion star. We exclude them from the primary sample because the unresolved binaries in hierarchical triplets can bias metallicity estimates by 0.1 dex (El-Badry et al. 2018). Excluding these pairs leaves us with 31 pairs of comoving stars.

We reduce the raw MIKE data with CarPy (Kelson 2003) which outputs wavelength-calibrated multi-order echelle spectra for the science targets and flats. After flat fielding, the spectra were normalized by fitting a cubic spline as a pseudo continuum for each order. We determine the pseudo continuum by iterative sigma clipping. We discard 100 pixels on both ends of each order because the blaze function is less well defined due to the lack of signal. Orders are then combined using a flux error weighted average. The sigma clipping process was tuned by comparing observations of four Gaia Benchmark stars to the online library (Blanco-Cuaresma et al. 2014b). We employ iSpec (Blanco-Cuaresma et al. 2014a) to determine the RV correction and remove cosmic rays.

The reduced spectra for two comoving pairs in the wavelength regime of 5384 ≤ λ ≤ 5400 Å are shown in Figure 2. In the top panel, we compare the spectra of both components for the two comoving pairs. The bottom panel shows the difference in spectra between the comoving components. Since we only target stellar twins, the difference predominantly comes from the difference in elemental abundances. The two spectra exhibit only minimal differences when they have similar chemistry. The careful selection of stellar twins allows us to derive elemental abundances with high fidelity for this homogeneous sample of comoving pairs.

Zoom In Zoom Out Reset image size

2.1. Positions and Velocities

We calculate the positions and velocities for each star using Gaia eDR3 astrometry and spectroscopic RVs from our observations. Note that we made the target selection with DR2 before eDR3 became available. However, since Gaia eDR3 improves Gaia DR2 astrometry, we decide to use eDR3 astrometry for the final determination and plots throughout this study. This update causes some objects within our main sample to appear outside the Δv_3D < 2 km s⁻¹ cut, even though our our main sample are selected to be within this boundary with the DR2 values. Table 2 summarizes the 3D seperation and Δv_3D calculations.

We will consider two approaches to calculate the 3D separations between members of comoving pairs. The first (method 1) uses the 6D phase space information to calculate the 3D separation directly using a Monte Carlo approach. We construct a multivariate normal distribution for each star with mean of the reported values and a covariance matrix using the 5D measured correlations in astrometry from Gaia eDR3 as well as the RV errors from the iSpec RV fit. We sample from this distribution 10,000 times and apply a transformation from equatorial to Galactic Cartesian coordinates to each sample. Finally, we calculate the 3D separation and velocity difference of the samples in Cartesian coordinates, taking the mean and standard deviation across samples as the point estimate and uncertainty for these quantities.

This method does not invoke additional assumptions about the geometry or separation distribution of pairs. However, for comoving pairs with 3D separations, ≲10⁵ au, the uncertainty in 3D separation is often significantly inflated by parallax uncertainties, making it impossible to accurately measure 3D separations using method 1. Consider two stars with identical true parallax, ϖ, and parallax uncertainty, σ_ϖ . In the limit of σ_ϖ /ϖ ≪ 1, the typical difference in apparent line-of-sight distance between the two stars is ${\rm{\Delta }}{d}_{\mathrm{apparent}}\sim \tfrac{1000\,\mathrm{pc}}{\left(\varpi /\mathrm{mas}\right)}\times {\sigma }_{\varpi }/\varpi$. For typical pairs in our sample, with ϖ ∼ 10 mas and σ_ϖ ∼ 0.02 mas, this translates to an apparent distance difference of ∼4 × 10⁴ au, which represents the amount by which the distance differences of typical pairs will be inflated by parallax uncertainties. For the pairs in our sample with the smallest ϖ and largest σ_ϖ , this value is ∼3 × 10⁵ au. This intuition also prompted us to divide the close comoving pairs and the far comoving pairs with a 1 pc (or 2 × 10⁵ au) threshold.

It is expected that many of the closest pairs in the sample are gravitationally bound wide binaries with true 3D separations smaller than this distance resolution afforded by the Gaia parallaxes. Parallax uncertainties effectively stretch these pairs out along the line of sight, in an effect that is analogous to the fingers-of-God distortion in cosmological surveys, where galaxy pairs and clusters are stretched out by redshift uncertainties. For such pairs, the point estimate of the 3D separation from method 1 will dramatically overestimate the true 3D separation. The projected separation (which is not subject to inflation by parallax uncertainties) will be much smaller than the point estimate of the 3D separation.

Figure 3 compares the 2D and apparent 3D pair separations for all 33 comoving pairs. The 3D separations are calculated using method 1, as described above. We also show the results separately using astrometry from DR2 and eDR3. At large separations, where the parallax uncertainties are small compared to the parallax differences, the 2D and 3D separations agree within a factor of a few, and the difference between the separations calculated using DR2 and eDR3 astrometry are minor. By contrast, pairs with 2D separations less than 10⁵ au fall far below the one-to-one line, with projected separations that are much smaller than the calculated 3D separations. For a majority of these pairs, the apparent 3D separation shrinks between DR2 and eDR3 astrometry. This reflects the fact that parallax uncertainties are smaller in eDR3 data, causing the separation inflation to be less severe. Crossmatching with El-Badry et al. 2021, hereafter E21) shows that almost all of these pairs have a high probability of being gravitationally bound wide binaries, ¹¹ suggesting that their 3D separations are indeed overestimated due to the parallax errors.

Zoom In Zoom Out Reset image size

This suggests that the separations of the close comoving pairs are dominated by parallax errors and that even Gaia eDR3 cannot resolve their line-of-sight distance difference. Therefore, for close pairs, we instead derive constraints on the 3D separation based on the projected separation, the assumption of random viewing angles, and a prior on the 3D separation distribution of wide binaries. This is describe in Appendix A. We refer to this method of estimating the 3D separation, which only relies on the projected 2D separation, as method 2. For method 2, the point estimate and 1σ confidence interval on the 3D separation, r_3D are ${r}_{3{\rm{D}}}={1.12}_{-0.11}^{+0.49}\times {r}_{2{\rm{D}}}$, where r_2D is the projected separation. Motivation for this expression is given in Appendix A. Table 2 summaries the separation and velocity differences between comoving pairs.

With different separation estimators for close pairs (thought to be binaries) and far pairs (thought to be unbound), it is necessary to decide where to draw the line between the two regimes. For simplicity, we use method 2 on all pairs for which the 3D separations calculated using method 1 are less than 2 × 10⁵ au (dashed blue line in Figure 3). It is possible that a few pairs wider than this (i.e. with 10⁵–10⁶ au), which also have projected separations significantly smaller than their calculated 3D separations, are also binaries. Improved parallax uncertainties in future Gaia data releases will help determine whether they are bound. But we have checked that separating the bound and unbound pairs with a 10⁶ au threshold leaves most of our qualitative results unaltered.

As most studies usually operate at the projected space instead of the 3D space adopted in this study, when comparing with the literature values (Ramírez et al. 2019, and Hawkins et al. 2020b), we calculate the 3D separations and velocity differences of their targets with the same process as detailed above. We adopt the latest Gaia eDR3 astrometry and the RV from the literature. RV data was not provided in Ramírez et al. (2019). Hence, we used Gaia eDR3 RVs for the Ramírez et al. (2019) data set.

We use the Brussels Automatic Code for Characterizing High accUracy Spectra (BACCHUS; Masseron et al. 2016) to determine stellar parameters (i.e., T_eff, log g, [Fe/H], v_micro) and abundances. An overview of the software is provided below. For a more comprehensive description we refer the reader to Section 2.2 of Hawkins et al. (2015); Masseron et al. (2016), and Appendix A.12 of Smiljanic et al. (2014).

BACCHUS fits observations using spectral synthesis. These spectra were constructed from MARCS model atmospheres (Gustafsson et al. 2008) using TURBOSPECTRUM (Plez 2012) for radiative transfer. MARCS models are calculated in 1D LTE. If the surface gravity (log g) is ≥3.0 dex, plane-parallel models are used, and spherical models otherwise. The atmospheric composition uses solar abundance from Grevesse et al. (2007) scaled by metallicity for most elements. MARCS models use separate abundance estimates for C, N, and O (see Gustafsson et al. 2008, Section 4 for more information). ¹² We use Gaia-ESO line list version 5 (Heiter et al. 2021) for atomic transitions. The line list includes hyperfine structure splitting for Sc ɪ vɪ, Mn ɪ, Co ɪ, Cu ɪ, Ba ɪɪ, Eu ɪɪ, La ɪɪ, Pr ɪɪ, Nd ɪɪ, and Sm ɪɪ. We also include molecular data for CH (Masseron et al. 2014), C₂, CN, OH, MgH (T. Masseron, private communication), SiH (Kurucz 1992), TiO, FeH, and ZrO (B. Pelz private communication). Fe ionization-excitation equilibrium is used in BACCHUS to derive effective temperature (T_eff), surface gravity, iron abundance ([Fe/H]), and microturbulent velocity (v_micro). Other broadening sources (e.g., rotation, instrument resolution, macroturbulence) are modeled by a Gaussian convolution.

Stellar parameters and the convolution kernel size are derived in BACCHUS using the param module. We first optimize for the convolution. The we determine T_eff by requiring a null trend,within 1σ of the fit, in Fe i abundance versus excitation potential. We determine log g by enforcing the abundances from Fe i and Fe ii lines to agree within the line-by-line rms and v_micro by enforcing that there is no trend in iron abundance to reduced equivalent width (i.e., EW/λ). The metallicity is taken as the mean of the Fe i lines relative to the Sun. The parameter fitting uses 94 Fe i lines and 32 Fe ii lines, and the details can be found in the full version of Table 3. The reported errors in T_eff, log g, and v_micro are calculated through the usual propagation of error from the determination of EW of individual lines. We note that, currently, BACCHUS does not estimate the covariances between stellar parameters.

Table 3. A Portion of Our Line List Selection

Identifier	Element	λ	$\mathrm{log}{gf}$	Reference Key	χ	$\mathrm{log}({A}_{X})$
		(Å)	(dex)		(eV)	(dex)
CM00A	Na ɪ	5682.630	−0.706	GESMCHF	2.102	6.053
CM00A	Na ɪ	5688.200	−0.404	GESMCHF	2.104	6.291
CM00A	Na ɪ	6154.220	−1.547	GESMCHF	2.102	6.181
CM00A	Na ɪ	6160.740	−1.246	GESMCHF	2.104	6.196
CM00B	Na ɪ	5682.630	−0.706	GESMCHF	2.102	6.103
CM00B	Na ɪ	5688.200	−0.404	GESMCHF	2.104	6.311
CM00B	Na ɪ	6154.220	−1.547	GESMCHF	2.102	6.212
CM00B	Na ɪ	6160.740	−1.246	GESMCHF	2.104	6.233
CM01A	Na ɪ	4751.820	−2.078	GESMCHF	2.104	6.558
CM01A	Na ɪ	5682.630	−0.706	GESMCHF	2.102	6.528

Note. A full machine-readable version, including abundances for each star for each line, is available online. The lines used will vary between stars because of the quality checks. χ is the excitation potential, and log(A_χ ) is the absolute abundance (before subtracting the Solar abundances) derived for this line. The solar abundances adopted are from Grevesse et al. (2007), except where described otherwise in Section 3. Reference keys can be matched with the full citation in Table 7 to determine the origin of the log gf values used.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

After optimizing the stellar parameters, abundances for 24 elements are measured in BACCHUS with the abund module. We measure abundances for the following elements: Li, C, Na, Mg, Al, Si, Ca, Sc, Ti, V, Cr, Mn, Fe, Co, Ni, Cu, Zn, Sr, Y, Zr, Ba, La, Nd, and Eu. The abund module synthesizes spectra at different [X/H], fixing the atmospheric model parameters to the best fitting stellar parameters. This model is generated through interpolation of nearby MARCS atmospheres; we determine the abundance for each line by minimizing the unweighted χ² of the observed and synthesized spectrum, and we repeat this process for all selected absorption lines for that species. The line list selection is given in Table 3. We select lines based on visual inspection of the synthesis compared with observed spectra for three stars, CM07A, CM11A, and CM25A. We select these stars because they are representative for different atmospheric parameters and rotational broadening. On top of that, for every star in our sample, selected lines must pass an internal quality check to be considered usable for that object. The internal quality check corresponds to a decision tree regarding whether the abundance output is physically plausible and constrained by the trial solutions (see Section 2.2 of Hawkins et al. 2015, for more details on the internal quality check). As a consequence of the internal quality check, the lines with measured abundances may vary between pairs. But we emphasize that for each pair, we adopt the same set of lines for the differential analysis.

We study the abundance differences using both line-by-line (differential) and non-line-by-line (nLBL) global methods. Both methods provide similar conclusions; however, the differential method generally leads to better precision because it mitigates systematic errors in the spectral models. We will adopt the abundances derived with the differential method throughout this study. But for completeness, we include results from both approaches in Table 4. For clarity, we denote abundances for individual lines with a subscript i. For example, [X_i/H] indicates an abundance of species X for line i, Δ[X_i/H] is the difference of species X using the abundance of line i between two stars. Species without subscripts are the overall estimate of the abundance from all corresponding spectral features.

Table 4. Stellar Parameters and Elemental Abundances Derived from BACCHUS

Component A ID	Component B ID	Identifier	Large ΔRV	${{\rm{T}}}_{{\mathrm{eff}}_{{\rm{A}}}}$	${\sigma }_{{T}_{{\mathrm{eff}}_{{\rm{A}}}}}$	${T}_{{\mathrm{eff}}_{{\rm{B}}}}$	${\sigma }_{{T}_{{\mathrm{eff}}_{{\rm{B}}}}}$	$\mathrm{log}{g}_{{\rm{A}}}$	${\sigma }_{\mathrm{log}{g}_{{\rm{A}}}}$	$\mathrm{log}{g}_{{\rm{B}}}$	${\sigma }_{\mathrm{log}{g}_{{\rm{B}}}}$	ξA
				(K)	(K)	(K)	(K)	(dex)	(dex)	(dex)	(dex)	(km s−1)
6160638833832941312	3471119180522314624	CM00	False	5809	9	5934	65	4.53	0.30	4.65	0.41	1.22
3613454637229633664	3613454637229634176	CM01	False	6058	25	6002	56	4.29	0.16	4.14	0.28	1.30
3639520621950395904	3639520621950395776	CM02	False	5962	11	5966	11	4.38	0.22	4.41	0.23	1.06
5897704951769034368	5897785667103983616	CM03	False	6036	37	6121	28	4.34	0.22	4.28	0.21	1.31
6224633983987510528	6224633983987511552	CM04	False	5763	26	5822	11	4.51	0.14	4.50	0.20	0.90

${\sigma }_{{\xi }_{{\rm{A}}}}$	ξB	${\sigma }_{{\xi }_{{\rm{B}}}}$	Δ[Fe/H]	${\sigma }_{{\rm{\Delta }}{[\mathrm{Fe}/{\rm{H}}]}_{\mathrm{Total}}}$	Δ[Na/H]	${\sigma }_{{\rm{\Delta }}{[\mathrm{Na}/{\rm{H}}]}_{\mathrm{Total}}}$	Δ[Mg/H]	${\sigma }_{{\rm{\Delta }}{[\mathrm{Mg}/{\rm{H}}]}_{\mathrm{Total}}}$	Δ[Al/H]	${\sigma }_{{\rm{\Delta }}{[\mathrm{Al}/{\rm{H}}]}_{\mathrm{Total}}}$	Δ[Si/H]
(km s−1)	(km s−1)	(km s−1)	(dex)	(dex)	(dex)	(dex)	(dex)	(dex)	(dex)	(dex)	(dex)
0.04	1.31	0.05	0.04	0.02	0.03	0.02	0.02	0.02	−0.02	0.02	0.01
0.04	1.33	0.04	0.02	0.02	0.05	0.02	0.03	0.02	0.03	0.01	0.01
0.04	1.05	0.04	0.00	0.02	0.00	0.02	0.01	0.02	−0.02	0.01	0.01
0.04	1.34	0.04	−0.08	0.02	−0.11	0.02	−0.06	0.02	−0.12	0.02	−0.08
0.04	0.99	0.04	−0.01	0.02	−0.00	0.03	0.01	0.01	0.00	0.01	−0.03

${\sigma }_{{\rm{\Delta }}{[\mathrm{Si}/{\rm{H}}]}_{\mathrm{Total}}}$	Δ[Ca/H]	${\sigma }_{{\rm{\Delta }}{[\mathrm{Ca}/{\rm{H}}]}_{\mathrm{Total}}}$	Δ[Sc/H]	${\sigma }_{{\rm{\Delta }}{[\mathrm{Sc}/{\rm{H}}]}_{\mathrm{Total}}}$	Δ[Ti/H]	${\sigma }_{{\rm{\Delta }}{[\mathrm{Ti}/{\rm{H}}]}_{\mathrm{Total}}}$	Δ[V/H]	${\sigma }_{{\rm{\Delta }}{[{\rm{V}}/{\rm{H}}]}_{\mathrm{Total}}}$	Δ[Cr/H]	${\sigma }_{{\rm{\Delta }}{[\mathrm{Cr}/{\rm{H}}]}_{\mathrm{Total}}}$	⋯
(dex)	(dex)	(dex)	(dex)	(dex)	(dex)	(dex)	(dex)	(dex)	(dex)	(dex)	⋯
0.02	−0.01	0.05	0.04	0.06	0.07	0.03	0.09	0.04	0.05	0.03	⋯
0.02	0.01	0.03	0.02	0.05	0.02	0.03	0.01	0.03	0.06	0.03	⋯
0.02	−0.02	0.04	0.01	0.06	−0.00	0.03	−0.01	0.04	0.01	0.03	⋯
0.02	−0.06	0.03	−0.10	0.05	−0.07	0.03	−0.09	0.04	−0.09	0.03	⋯
0.03	−0.00	0.05	−0.04	0.07	0.00	0.03	0.05	0.04	0.03	0.04	⋯

Note. A and B subscripts denote the two components from the comoving pairs. ξ denotes the microturbulence. The electronic table contains measurements from nLBL and differential methods. All the differential measurements (e.g., Δ[Fe/H], σ_Δ[Fe/H]) are measured with differential method. But we also provide individual estimates of each component with the nLBL method, indicated with the subscripts A and B. By definition, the differential analysis only applies to the differences of the two components, and there will be no individual values with the differential method. The full table also separately lists the errors from differential stellar parameters (denoted with a subscript θ) and the statistical line-to-line scatter (denoted with a subscript line), but below we only show the total errors (adding in quadrature). ΔRV refers to the difference in RV measured from our observations and Gaia eDR3. The large ΔRV column indicates whether the system was flagged for potentially being a hierarchical triplet with an unresolved binary.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

In the differential approach, we compute Δ[X_i/H] for each transition observed in both stars. The abundance difference is taken as the median of these values (i.e., Δ[X/H] =median(Δ[X_i/H])). The median is used instead of the mean because it is more robust against outliers. As for the certainty estimation, we assume the statistical errors for Δ[X/H], i.e., std(Δ[X_i/H])/$\sqrt{\mathrm{number}\ \mathrm{of}\ \mathrm{lines}}$, as the standard error of the mean of the line abundance differences. However, if a single line remains after quality cuts, we take the error as 0.1 dex. This is 1/3 the step size used between trial values in [X/H] when calculating the abundance for each line.

Abundance errors from the uncertainty in T_eff, log g, and v_micro are propagated according to Hawkins et al. (2020b). For each parameter, we perturb the best-fit model and derive abundances for this perturbed model atmosphere. The difference between the abundances from the best-fit and perturbed model is the error introduced from that parameter. These abundance errors are added in quadrature with the line-by-line statistical errors for Δ[X/H] to determine the total error for an abundance measurement. For simplicity, and due to the limitations of the BACCHUS code, we choose to omit the covariances following Hawkins et al. (2020b).

We perform the error analysis on a subset of 30 representative stars distributed evenly over the stellar parameters of our data set. The offsets for T_eff, log g, and v_micro , where ±20 K, ±0.1 dex, and ±0.03 km s⁻¹, respectively. The values for T_eff and v_micro are the median 1σerrors reported by BACCHUS for our sample. We assume a Δ log g = 0.1 dex instead of the median reported errors from BACCHUS (0.20 dex) because we find that BACCHUS tends to overestimate the errors of log g due to some outlying Fe ii lines. We have separately analyzed these spectra with the Bayesian stellar parameter code LoneStar (Nelson et al in prep), which agrees with BACCHUS estimates on all quantities besides the log g errors, which are found to be ∼0.05 dex. We chose log g errors of 0.10 dex to be conservative since BACCHUS favors larger log g errors. We note however that our comparison between the close comoving pairs, far comoving pairs, and random pairs is robust and independent of our error estimates.

When we estimate the errors for individual observations, we adopt interpolation by inverse distance weighting (IDW) to estimate the abundance errors resulting from the uncertainties in the stellar parameters for the remainder of the data set. IDW interpolation approximates points outside the reference set by constructing a weighted average of points from the reference sets. These weights are the inverse rms distance. During the interpolation, we transform the stellar parameters onto the interval (−1, 1) to weigh each equally. We tested the accuracy of the interpolation using a leave-one-out cross validation for each element. The average error introduced was typical <2% of the true errors calculated for that particular star. Sr and Eu were the only two elements with larger relative errors at ≲5%. Nonetheless, these errors are at least one order of magnitude below the total error budget for each element and are therefore negligible.

Table 4 summarizes the derived stellar parameters and abundance measurements and uncertainties in this study. In addition to the total abundance errors, we also report the errors from the individual errors contributed by the stellar parameter error propagation and the line-to-line scatter of Δ[X/H], respectively. The stellar parameters are the main source of uncertainty. For 23 elements, the stellar parameters account for >90% of total error. For Sr, the stellar parameters account for ∼70% of the total errors. Of the remaining species, Li and Eu typically only have a single line measured and are therefore subjected to the 0.1 dex error floor on the abundance scatter. Only in these rare cases, the line-to-line error is more dominant than the stellar parameters uncertainty. Finally, for completeness, although not shown, the nLBL abundances and errors are also provided in the full machine-readable table.

As discussed in the introduction, this papers main goal is to assess whether or not the comoving pairs in this study, especially for the far comoving pairs, are conatal systems. Since conatalitly cannot be assessed directly outside star-forming regions, we resort to studying the chemical homogeneity (i.e., Δ[X/H]). Figure 4 displays the difference in metallicity between comoving pairs (blue) as a function of pair separation. Although not shown, similar trends are seen for the majority of the elements studied.

Zoom In Zoom Out Reset image size

For comparison, we also include the wide binaries results from Ramírez et al. (2019, and Hawkins et al. (2020b). We note that Ramírez et al. (2019) includes results for 12 systems, 11 of which are from other studies (i.e., Ramírez et al. 2010; Liu et al. 2014; Mack & Schuler 2014; Ramírez et al. 2014; Tucci Maia et al. 2014; Teske et al. 2016; Saffe et al. 2016; Mack & Stassun 2016; Saffe et al. 2017; Reggiani & Meléndez 2018; Oh et al. 2018). In the remainder of the paper, we will exclude the system from Mack & Schuler (2014) from comparisons because it appears to be a hierarchical triplet with an unresolved binary, and use Ramírez et al. 2019 as a shorthand for the remaining 11 systems. We note that since the sample in Ramírez et al. (2019) is a heterogeneous sample, the comparison with Hawkins et al. (2020b) might be more relevant. Furthermore, the samples in Hawkins et al. (2020b) have spectra with spectra and were analyzed in the same way as this study. We also compare our data to Andrews et al. (2019); however, they did not provide individual measurements, so our comparison with this work is limited to summary statistics.

As previously discussed, we split our comoving sample into the close (i.e., primarily wide binaries) and far comoving pairs to compare with these works. Close comoving pairs are defined as those having separations below 1 pc, i.e. 2 × 10⁵ au, (with 17 pairs in this study), and far comoving pairs have separations above 1 pc (14 pairs). For the close pairs, we use the geometric prior when calculating the separations (see Section 2.1) because the parallax error dominates. Our close comoving pairs sample (i.e., wide binaries) show a similar degree of chemical homogeneity as other studies on wide binaries with a dispersion of 0.05 dex in Δ[Fe/H], which is expected. Compared to the close comoving pairs, we find that the far comoving pairs are less chemically homogeneous, with a Δ[Fe/H] dispersion of 0.08 dex; however, the far comoving pairs are substantially more homogeneous than the random pairs (0.23 dex), which we will further explore in the next section. Table 5 summarizes the scatter in Δ[Fe/H] comparing our data with prior works.

4.1. Chemical Homogeneity of Comoving Pairs

To contextualize how homogeneous our comoving population is, we simulate pairs of random field stars. In a similar fashion to Hawkins et al. (2020b), we create a set of random pairs by assigning every star to another star, which is not its comoving partner. Since the comoving stars are selected to be stellar twins, for better comparison, we create random pairs, which minimize the difference in T_eff and log g (labeled T_eff/log g pairs) to mitigate systematic effects. Since T_eff and log g have different units, when minimizing the difference, we divide both quantities by the standard deviation of the pairwise difference for all non-comoving pairs of stars in our sample. The random pairs are shown in Figure 5. Across all elements, the average dispersion in Δ[X/H] for the random pairs is 0.24 dex, significantly larger than the dispersion of both the close comoving pairs (0.05 dex) and the far comoving pairs (0.10 dex).

Zoom In Zoom Out Reset image size

To better illustrate the difference in chemical homogeneity between these groups, similar to Andrews et al. (2019), we visualize the homogeneity of our sample using ∣Δ[X/H]∣/σ_Δ[X/H], which will be abbreviated as Δ/σ. The Δ here is the difference in elemental abundances between the two components, and the σ is the quadrature sum of the uncertainties from the two components. Figure 6 shows the cumulative distribution function (CDF) of Δ[X/H] normalized by the measurement error for the comoving pairs (violet/solid and dashed) and random pairs (teal/solid). For reference, each panel includes the CDF of Δ/σ drawn from the unit Gaussian, i.e., what a chemical homogeneous would look like at the level of our current measurement uncertainties (assuming that our uncertainty estimations are accurate).

Zoom In Zoom Out Reset image size

There is a clear difference in the random pairs distribution and comoving pairs, both close and far, largely attributed to the larger chemical dispersion seen in the random pairs. The difference in the CDF between the comoving and random pairs also varies between elements. This variation is primarily caused by the difference in our measurement precision for those elements. On the one hand, for better-measured elements, such as Fe, Mg, Al, and Si, the comoving pairs exhibit a more distinct chemical homogeneity than the random pairs. On the other hand, unsurprisingly, for less well-measured elements, the difference in chemical homogeneity is less distinct, as shown in the CDFs for La, Nd, and Eu. We do not include Li in this plot because its surface abundance has a strong dependence on effective temperature and is therefore not a useful indicator for conatality (see Gray 2005; Ramírez et al. 2012; Hawkins et al. 2020b).

For most elements, the close comoving pairs are consistent with a chemically homogeneous distribution at our measurement precision, with the CDF for 16 elements exceeding 96% by Δ/σ = 4. In contrast, the random pairs are less homogeneous, with the CDF for 13 elements ≲50% at Δ/σ = 4.

As already seen in Figure 4, the far comoving pairs are more heterogeneous than the close comoving pairs for most elements, albeit still being more homogeneous than the random pairs. The slightly larger chemical inhomogeneity is unsurprising because we expect some of the far comoving pairs might be contaminated by chance alignment pairs (e.g., Kamdar et al. 2019a), which we will quantify in the next section.

In this study, we focus on the abundance [X/H] instead of the abundance ratio [X/Fe]. However, all elements correlate with iron in [X/H]. Therefore, the difference between random and close comoving pairs in Δ/σ for [X/H] may only reflect a metallicity difference. While [X/H] (or metallicity alone) are telltale signs of conatality (or at least coeval), studying the abundance ratio [X/Fe] might, in principle, constitute a more stringent test and might reveal subtle physics. Although not shown, we tried to examine Δ/σ for [X/Fe] ¹³ . We find little difference between the comoving and random pairs of stars in the [X/Fe] space. This is not unexpected because the residual variance in [X/Fe] is much more subtle than in [X/H]. For example, Ting & Weinberg 2021 show that when one subtracts the mean chemical track, the residual correlation in [X/Fe] is small, with a signal ≲0.01–0.02 dex. Therefore, to see such a signal, we would need to measure [X/Fe] better than ∼0.01 dex. With the current pipeline adopted in this study, we found that achieving 0.01 dex remains challenging. Therefore, we only focus on [X/H]. Future studies with even higher quality spectra and/or better analysis will be needed to shed light on these subtle variations.

4.2. Conatal Fraction for Separations >2 × 10⁵ au

The close comoving pairs (wide binaries) at separations below 2 × 10⁵ au are commonly believed to be conatal as chance capture is unlikely. However, beyond such separation, the conatality for the far comoving pairs remains unexplored. Armed with our high-resolution spectroscopic data from these pairs, in this section, we will constrain how conatal the far comoving pairs with separations >2 × 10⁵ au are relative to the close comoving pairs. Showing the conatality of far comoving pairs can have far-reaching consequences; if these comoving pairs are chemically homogeneous and conatal, this would significantly expand the sample of calibrators beyond open clusters and wide binaries.

As we have seen in the previous section, these far comoving pairs are more heterogeneous than the close comoving pairs/wide binaries but more chemically homogeneous than the random pairs. One explanation for the slightly larger chemical inhomogeneity is that the far comoving pairs are a mixture of conatal pairs and chance alignments. The interlopers create a tails of large Δ [X/H]. If that is the case, we should be able to constrain the conatality fraction through these far comoving pairs chemical homogeneity, which is what we will attempt next.

We only model Δ [Fe/H] because iron has the most precise measurements of any element in our sample, and as we have argued, most elemental abundances only trace [Fe/H] at the current precision. The basic idea is that if the Δ [Fe/H] of the far comoving pairs are contributed by the two populations, a conatal population that resembles the wide binaries population and the random pairs interlopers, we would expect the distribution of Δ [Fe/H] would also be a mixture of the two. More specifically, we model the probability density function (PDF) of ∣Δ[Fe/H]∣/σ_Δ[Fe/H] (or in short, Δ/σ, in the following) for the far comoving pairs as a mixture model using the weighted average of the close comoving pairs' Δ/σ PDF and the random pairs' Δ/σ PDF. The best-estimated weight signifies the fraction of pairs consistent with the close comoving pairs among the far comoving populations, hence its conatal fraction. Due to our small sample size (∼10 pairs), we expect considerable uncertainty for estimating the conatal fraction because of the sampling noise. Therefore, to properly model the sampling noise, we also estimate the uncertainty of the estimated conatal fraction by bootstrapping the sample.

Operatively, first, we sample, with replacement, the data distribution of Δ/σ for the close, far, and random pairs. The distributions of points drawn from the close comoving pairs and random pairs are then separately approximated as normal distributions through maximum likelihood estimation (MLE). Subsequently, we fit for the weight of the distribution of the far comoving pairs, also through MLE, assuming that the Δ/σ distribution of the far comoving pairs is a weighted mixture of the ones determined for the close moving pairs and the random pairs. The bootstrapping process is repeated 1000 times. The bootstrapping produces a distribution of weights. The median and standard deviation are taken as the best estimate and uncertainty of the conatal fraction, respectively.

In Figure 7, we demonstrate the nominal fit without the replacement. For this estimate, we only consider comoving pairs with separations 2 × 10⁵–10⁷ au, excluding the three pairs with separations >10⁷ au, leaving us with 11 pairs (out of 14 pairs). As we will see in the following, these three pairs are clearly interlopers (based on predictions from simulations). The Δ/σ CDF of the close moving pairs is shown in violet, the far comoving pairs in violet/black, and the random pairs in teal. A continuum of mixture model CDFs is displayed in the background as a blue/yellow gradient ranging from close:random ratios of 1:4 to 9:1, or equivalently, a conatal fraction 20%–90%. Besides, we highlight some fractions with white dotted contours to guide the reader. We find the conatal fraction for the 11 far comoving pairs with separations 2 × 10⁵–10⁷ au to be 73 ± 22%. In short, our result implies that ∼75% of the far comoving pairs are conatal based on their metallicity measurements.

Zoom In Zoom Out Reset image size

Conatal stars are the central hallmarks in modern-day galactic archaeology. Absolute standards are hard to come by (Jofré et al. 2014); Blanco-Cuaresma et al. 2014b); Heiter et al. 2015), but if we know the two stars are conatal, the two components serve as each other references. This is why open clusters and wide binaries have always been the gold standard in calibrators to refine stellar models, to study exoplanets, gas mixing, and to calibrate surveys. More recently, the high precision astrometry measurements from Gaia now allow us to identify wide binaries at high fidelity, which further propels the study of wide binaries. However, most studies of wide binaries (e.g., El-Badry et al. 2021) often assume a dichotomy between wide binaries and chance alignments, and this picture is clearly too simplistic. Simulations (Kamdar et al. 2019a, 2019b) have suggested comoving stars from disrupted star clusters could have a high probability of being conatal out to separations of ${ \mathcal O }({10}^{6})$ au. If this is correct, unbound comoving stars could be used in similar ways as wide binaries and open clusters. Based on modeling of Gaia data by Y.-S. Ting et al (2020, in preparation), we expect to find roughly 3300 conatal pairs with 3D separations ≲2 × 10⁶ au (10 pc) and Δv_3D < 2 km s⁻¹.

In this study, we perform the first homogeneous study of comoving pairs that span five orders of magnitude in separation. A homogeneous sample allows us to study the chemical homogeneity of the comoving stars with different separations consistently. We found that the close comoving pairs with separations <2 × 10⁵ au have a Δ[Fe/H] scatter of 0.05 dex. These values are comparable to the typical values seen in open clusters and other wide binaries studies. In particular, we opt to perform the same analysis as Hawkins et al. (2020b). With similar resolution and S/N, Hawkins et al. (2020b) showed a dispersion of 0.05 dex for wide binaries, largely consistent with this study. On top of that, we show that the abundance differences for almost all elements are consistent with the measurement uncertainties, further validating that close comoving pairs/wide binaries are conatal. We note that some of the Δ[Fe/H] scatter seen in the systems from Ramírez et al. (2019) could be a consequence of a selection effect. Many of those systems were published in the context of planet engulfment signatures. Consequently, the systems that show the engulfment signals (e.g., Kronos and Krios Oh et al. 2018) will also have larger metal differences. Since exoplanet accretion changes the surface abundances, these systems may not be representative of the true chemical homogeneity of conatal systems.

More importantly, our study expanded on previous works by also examining comoving pairs with separation >10⁶ au. The key result of our study is that comoving pairs at these larger separations are still significantly more chemical homogeneous than the random pairs, albeit with a slightly larger dispersion in chemistry, with a scatter Δ[Fe/H] of 0.09 dex. The random pairs have a scatter of 0.23 dex. If we assume that the far comoving pairs are comprised of both the conatal stars and chance alignments, we estimate that about 73 ± 22% of pairs with separations 2 × 10⁵–10⁷ au are conatal.

A caveat with this estimation is that, in this study, we select stellar twins with a median ∣ΔT_eff∣ of 105 K. The selection of stellar twins enables the mitigation of any potential systematics in terms of the spectral modeling by performing a line-by-line differential study. Conveniently, from their construction, the simulated random pairs have a smaller difference in stellar parameters (with a median ∣ΔT_eff∣ of 47 K), so we would expect these systematic effects to be similar or weaker in the random population. However, for completeness, we investigated whether or not the more significant dispersion in metallicity difference is real or simply due to systematic errors when deriving abundances from stars that are more different in stellar parameters.

To dissect that, we looked into any potential biases in Δ[Fe/H] as a function of the difference in stellar parameters. We found weak trend between Δ[Fe/H] and ${\rm{\Delta }}{{\rm{T}}}_{\mathrm{eff}},{\rm{\Delta }}\mathrm{log}g,$ and Δv_micro. For example, Δ[Fe/H] shows a positive correlation with ΔT_eff, with a gradient of 2.3 × 10⁻⁴dex K⁻¹. We attempted to account for that by fitting a multivariate linear regression between Δ[Fe/H] as a function of ${\rm{\Delta }}{T}_{\mathrm{eff}},{\rm{\Delta }}\mathrm{log}g,$ and Δv_micro. We found that this process typically leads to a correction of 0.002 dex in Δ[Fe/H], which is negligible for our study. We conclude that it is unlikely the differences in chemistry between the comoving and random pairs are attributable to this correlation. We opted not to remove these correlations because it is hard to determine the exact causal direction of this correlation.

Recall that the observations in this study are largely motivated by the Kamdar et al. (2019a) simulation, in which the authors argued that far comoving pairs are also conatal. So, how does our study compare with the simulation? In Figure 8, we compare our data with the Kamdar et al. (2019a) simulation. Our comoving pairs are represented as red circles, and the symbol size shows the metallicity difference of the two stars for individual pairs. The black dashed lines indicate the escape velocity of bound systems, assuming equal-mass companions, and with the total masses of 2 and 5M_⊙ respectively. The simulation does not extend below 10⁵ au because the simulations did not attempt to model wide binaries (only disrupted star clusters).

Zoom In Zoom Out Reset image size

Our results exhibit excellent agreement with the simulations; our sample shows a high degree of chemical homogeneity in the regime where the simulations predict to have a high conatality rate (the region in yellow). In the phase space regime where chance alignments should dominate (in blue background), our observations also show the largest metallicity differences (largest symbol sizes). Also more quantitatively, the simulation predicts that ∼80% of comoving pairs are conatal provided Δv_3D < 2 km s⁻¹ and the pair separation is between 2 × 10⁵ and 10⁷ au. Our conatality fraction estimate of 73 ± 22% for separations is well aligned with the simulation's predictions.

In addition, besides our (31 pairs of the main sample), as discussed in the target selection, we also observed two pairs with large Δv_3D (>2 km s⁻¹) as the control sample. One of the two pairs has a spatial separation of 2 × 10² au. Even though they are not within Δv_3D < 2 km s⁻¹, this pair is clearly a bound binary, and therefore is conatal. The other pair has a large separation (2 × 10⁶ au) and a large velocity difference (4 km s⁻¹). The simulation predicted this pair is firmly a chance alignment pair as disrupted star cluster members are unlikely to exhibit such a large velocity difference. As shown in Figure. 8, our abundance measurement indeed concurs with this picture, the former has a metallicity difference that is typical of that of wide binaries, and the latter has a dispersion of 0.24 dex, consistent with what we expect from a random pair. With this in mind, we caution that all the results presented in this study, including the conatal fraction, only strictly applies to stars with Δv_3D < 2 km s⁻¹.

While the agreement is encouraging, due to the small sample size, unfortunately, we could only settle on two separation bins—those with separations <2 × 10⁵ au and those with 2 × 10⁵–10⁷ au. It is clear from the simulation that there is no sharp transition between the conatal to random pairs. But instead, the regime of purely conatal pairs transitions gradually toward the regime dominated by chance alignments. The transition depends critically on many Galactic properties, including how stars form and disperse and the cluster mass function. Also, as explored in Kamdar et al. (2019a), the number density of comoving stars provide telltale signs of the number density of Galactic perturbers, such as giant molecular clouds in the Milky Way disk, as they disrupt comoving pairs. In short, understanding this transition with critical and refining this boundary region with a larger sample will no doubt offer insights into the formation and dispersion of star clusters and other substructures in the Milky Way.

The exact value of the division between close pairs (i.e., wide binaries) and far pairs (presumed unbound) is somewhat arbitrary. Bound pairs are not expected to exist at separations beyond 1–2 pc because at these separations the Galactic tidal field is stronger than the internal acceleration within a binary (e.g. Jiang & Tremaine 2010). At separations of order 1 pc, there will always be some ambiguity between still-bound and recently dissolved pairs, and the reliability with which bound pairs can be identified will depend on the precision of the astrometry (e.g. El-Badry et al. 2021). A possible alternative to adopting a strict threshold at S_eDR3 = 2 × 10⁵ au would be to only consider as binaries pairs that are members of a predetermined wide binary catalog, such as the one produced by El-Badry et al. (2021). The disadvantage of this approach is that such catalogs are generally not 100% complete. Increasingly precise astrometry from future Gaia data releases will make it possible to more reliably classify pairs with separations of order 1 pc. For now, we note that because of our small sample size, our estimates of the conatal fraction are considerable and are dominated by small number statistics. Small changes in the adopted boundary between close and far pairs thus lead to changes in our inferred conatal fraction that are comfortably within our reported uncertainties. For completeness, we provide a summary of the effects of choice of separations for the far pairs in Table 6.

Finally, in this study, we assume that the difference in chemical homogeneity for the far comoving pairs compared to the close conatal pairs stems from random pairs interlopers. However, another potential source for these differences in chemical homogeneity could be from the ISM being less mixed at the largest scale. At the ${ \mathcal O }({10}^{7})$ au scale, ISM can have abundance scatter of up to 0.3 dex (e.g. Sanders et al. 2012). If stars formed in a poorly mixed cloud, they might inherit this scatter. We deem such a scenario unlikely because hydrodynamic simulations have found that, in a single star-forming region, even small amounts of turbulence can cause a factor of 5 reductions in abundance to scatter in stars compared to the progenitor ISM (Feng & Krumholz 2014). And this chemical reduction scatter applies to scales as large as ${ \mathcal O }({10}^{6})\mbox{--}{ \mathcal O }({10}^{7})$ au, the maximum scale probed in this study.

That said, if the stars formed in a filamentary structure at different star-forming regions, the ISM might be less well mixed. For example, Hawkins et al. 2020a argue the Pisces-Eridanus stellar stream could have been produced from a stellar filament. They find a metallicity range of ∼0.2 dex for the stream, which is consistent with our metallicity difference in the far comoving pairs. If the chemical inhomogeneity is due to the ISM physics instead of random interlopers, our conatality fraction estimated would be a conservative limit. The far comoving pairs would have an even higher conatal fraction than what was is inferred here.

We obtained a homogeneous sample of high-resolution (R ≃ 40,000), high S/N (∼150 pixel⁻¹) observations of 62 FG-type stars residing in 31 comoving pairs with separations that span five orders of magnitude in separation, between 10² au and 10⁷ au. Previous studies have mostly been restricted to wide binaries with a separation <10⁶ au, and we investigate if the (unbound) comoving stars with separations = 2 × 10⁵–10⁷ au are conatal. We measured the stellar atmospheric parameters and chemical abundances for 24 species, including covering the different nucleosynthetic pathways. The derived elemental abundances: Li, C, Na, Mg, Al, Si, Ca, Sc, Ti, V, Cr, Mn, Fe, Co, Ni, Cu, Zn, Sr, Y, Zr, Ba, La, Nd, and Eu.

We separate our sample into the classical wide binaries regime (i.e., close comoving pairs with separations < 2 × 10⁵ au), the far comoving pairs (separations >2 × 10⁵ au), as well as the random field pairs created through random pairings of our observations. We find the wide binaries are significantly more chemically homogeneous than field stars in [X/H] (see Figure 6) with a typical scatter in Δ[Fe/H] of 0.05 dex. Our results agree with previous works on wide binaries and comoving pairs (see, e.g., Andrews et al. 2019; Hawkins et al. 2020b). For most elements in the close comoving pairs, the dispersion in [X/H] is consistent with the measurement errors, implying that the close comoving/wide binaries pairs are chemically homogeneous at our measurement precision.

We demonstrate that the far comoving pairs with separations of 2 × 10⁵–10⁷ au also exhibit substantial chemical homogeneity (0.09 dex in the [Fe/H] scatter) compared to random pairs (0.23 dex). Nonetheless, the far comoving pairs are less homogeneous than the close comoving/wide binary population. If we assume that the far comoving pairs comprise a mixture of the conatal and random chance alignment population, modeling the distribution of Δ[Fe/H] as a mixture of these two populations implies that about 73 ± 22% of the unbound comoving stars with separations 2 × 10⁵–10⁷ au and Δv_3D < 2 km s⁻¹ are conatal. This conatality fraction is in excellent agreement with the predictions from simulations from Kamdar et al. (2019a).

Our study implies that most comoving stars are conatal, even though they are well separated. As well-separated pairs of stars are more common than wide binaries, this new population of clusters of two enables many windows for studies that we have thus far restricted to the wide binaries and open clusters. Harnessing these vastly abundant comoving pairs will have broad applications, ranging from calibrating surveys to understanding star formation, planet engulfments, and beyond.

T.N. and K.H. have been partially supported by a TDA/Scialog (2018–2020) grant funded by the Research Corporation and a TDA/Scialog grant (2019–2021) funded by the Heising-Simons Foundation. T.N. and K.H. acknowledge support from the National Science Foundation grant No. AST-1907417. Y.S.T. is grateful to be supported by the NASA Hubble Fellowship grant No. HST-HF2-51425.001 awarded by the Space Telescope Science Institute. K.H. is partially supported through the Wootton Center for Astrophysical Plasma Properties funded under the United States Department of Energy collaborative agreement No. DE-NA0003843. A.P.J. acknowledges support from a Carnegie Fellowship and the Thacher Research Award in Astronomy. H.K. acknowledges support from the DOE CSGF under grant No. DE-FG02-97ER25308. The authors thank Carnegie Observatory for granting us the observing time to conduct this study.

This research made use of the SIMBAD database, operated at CDS, Strasbourg, France (Wenger et al. 2000), and NASAs Astrophysics Data System Bibliographic Services. This work also made use of data from the European Space Agency mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC,https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

Facilities: Magellan/Clay Telescope - , Simbad - .

Software: astropy (Astropy Collaboration et al. 2013, 2018), NumPy (Harris et al. 2020), iPython (Perez & Granger 2007), Matplotlib (Hunter 2007), Galpy (Bovy 2015), SciPy (Virtanen et al. 2020), Photutils (Bradley et al. 2020), BACCHUS (Masseron et al. 2016), CarPy (Kelson 2003), topcat (Taylor 2005), iSpec (Blanco-Cuaresma et al. 2014a).

Consider a pair of stars with true 3D separation r_3D. If the pair is viewed along a random viewing angle, the probability distribution of its projected separation r_2D can be calculated via a straightforward geometric argument (e.g. Nottale & Chamaraux2018):

$$\begin{eqnarray}&&p\left({r}_{2{\rm{D}}}| {r}_{3{\rm{D}}}\right)=\displaystyle \frac{1}{{r}_{3{\rm{D}}}}\displaystyle \frac{{r}_{2{\rm{D}}}/{r}_{3{\rm{D}}}}{\sqrt{1-{\left({r}_{2{\rm{D}}}/{r}_{3{\rm{D}}}\right)}^{2}}}.\end{eqnarray} \tag{ A1 }$$

For the pairs in our sample suspected to be binaries, we have a precise measurement of r_2D and wish to constrain r_3D. By Bayes' rule, the probability distribution of r_3D given a measurement or r_2D depends both on the measured r_2D, and on the prior on r_3D:

$$\begin{eqnarray}&&p\left({r}_{3{\rm{D}}}| {r}_{2{\rm{D}}}\right)\propto p\left({r}_{2{\rm{D}}}| {r}_{3{\rm{D}}}\right)p\left({r}_{3{\rm{D}}}\right).\end{eqnarray} \tag{ A2 }$$

The separation distribution of wide binaries is well-constrained observationally to be $p\left({r}_{3{\rm{D}}}\right)\sim {r}_{3{\rm{D}}}^{-3/2}$ at r_3D > 500 au (e.g. Andrews et al. 2017; El-Badry & Rix 2018), so we adopt this as the prior for the 3D separation distribution of pairs thought to be binaries. In this case, Equation A2 becomes

$$\begin{eqnarray}&&\begin{array}{l}p\left({r}_{3{\rm{D}}}| {r}_{2{\rm{D}}}\right)\\ \ =\ \left\{\begin{array}{ll}\displaystyle \frac{2}{\sqrt{\pi }}\displaystyle \frac{{\rm{\Gamma }}\left(7/4\right)}{{\rm{\Gamma }}\left(5/4\right)}\displaystyle \frac{1}{{r}_{3{\rm{D}}}}\displaystyle \frac{{\left({r}_{2{\rm{D}}}/{r}_{3{\rm{D}}}\right)}^{5/2}}{\sqrt{1-{\left({r}_{2{\rm{D}}}/{r}_{3{\rm{D}}}\right)}^{2}}}, & {r}_{3{\rm{D}}}\gt {r}_{2{\rm{D}}}\\ 0 & \mathrm{otherwise}\end{array},\right.\end{array}\end{eqnarray} \tag{ A3 }$$

where Γ represents the Gamma function, and the coefficient is derived from the normalization condition. The red line in Figure 9 shows this distribution, with r_2D = 10⁴ au adopted for concreteness. ¹⁴ For this distribution, the median value of r_3D/r_2D is 1.12, and the middle 68.2% range is (1.01, 1.61). Given a measurement of r_2D, we thus report ${r}_{3{\rm{D}}}={1.12}_{-0.11}^{+0.49}\times {r}_{2{\rm{D}}}$ if the pair is suspected to be a binary.

Zoom In Zoom Out Reset image size

Of course, other choices are possible for the prior p(r_3D) when the separation exceeds what we expect from the binary population. The black line in Figure 9 corresponds to a flat prior, $p({r}_{3{\rm{D}}})=\mathrm{const}$. Such a separation distribution might be expected for conatal pairs that are no longer bound, if the birth and dissolution rates of such pairs are constant. In this case, the median and 1σ range in r_3D would be ${r}_{3{\rm{D}}}={1.41}_{-0.38}^{+2.63}\times {r}_{2{\rm{D}}}$, still implying that most pairs have 3D separations within a factor of a few of their 2D projected separation.

Finally, the dotted cyan line in Figure 9 shows the results of assuming $p({r}_{3{\rm{D}}})\propto {r}_{3{\rm{D}}}^{2}$. This is the distribution expected for random pairings of isotropically distributed stars. In this case, the probability distribution is normalizable only if we assume that there are no pairs beyond a 3D separation ${r}_{\max }$. The conditional probability distribution in this case is

$$\begin{eqnarray}p\left({r}_{3{\rm{D}}}| {r}_{2{\rm{D}}}\right)=\left\{\begin{array}{ll}\displaystyle \frac{1}{\sqrt{\left(1-{r}_{2{\rm{D}}}^{2}/{r}_{3{\rm{D}}}^{2}\right)\left({r}_{\max }^{2}-{r}_{2{\rm{D}}}^{2}\right)}}, & {r}_{3{\rm{D}}}\gt {r}_{2{\rm{D}}}\\ 0 & \mathrm{otherwise}\end{array},\right.\end{eqnarray} \tag{ A4 }$$

which approaches a constant value at r_3D ≫ r_2D. Intuitively, what happens in this case is that the larger number of possible pairs at large r_3D (which scales as ${r}_{3{\rm{D}}}^{2}$) exactly compensates for the smaller number of sightlines along which a pair can be viewed to have a given projected separation (which scales as ${r}_{3{\rm{D}}}^{-2}$). In Figure 9, we adopt ${r}_{\max }={10}^{6}\,\mathrm{au}$.

To summarize, we calculate constraints on the 3D separation of pairs thought to be wide binaries based on their projected 2D separation using Equation A3. This yields a median and 1σ range of ${r}_{3{\rm{D}}}={1.12}_{-0.11}^{+0.49}\times {r}_{2{\rm{D}}}$, which at small r_2D is much more constraining than the constraint on r_3D that considers the distances to both stars independently (method 1). The primary assumption of this constraint is that the separation distribution of wide binaries falls off as $p({r}_{3{\rm{D}}})\propto {r}_{3{\rm{D}}}^{-3/2}$. It is not appropriate for pairs that are not binaries (unbound comoving pairs or chance alignments) because their separation distribution is expect to increase with separation due to the different priors. We thus only apply this correction to pairs with separation <2 × 10⁵ au.

A large part of this work is made possible by numerous heroic efforts by various groups who perform the pain-stacking tasks of curating/calibrating the atomic line list used in BACCHUS. Those efforts are unfortunately under appreciated in the literature. As such, albeit long, we decide to include the references to all original sources of the line list adopted in this study in Table 7.