Stars that Move Together Were Born Together

I. Introduction

Close stars moving together in the Galaxy may hold valuable clues related to star formation (Reipurth & Mikkola 2012; Parker et al. 2011, e.g.,) and the dynamical history of the Galaxy (Weinberg et al. 1987; Monroy-Rodríguez & Allen 2014, e.g.,). Pairs of stars that are close together ($\lesssim 1$ pc) are likely gravitationally bound (Jiang & Tremaine 2010). Pairs further apart ($\gtrsim 1$ pc) are less likely to be gravitationally bound and may be associated with dissolving clusters (Kouwenhoven et al. 2010), thereby allowing us to study cluster disruption and the star formation history of the Milky Way (Bland-Hawthorn et al. 2010, e.g.,).

The study of co-moving pairs using proper motions has a rich history (Poveda et al. 1994; Chanamé & Gould 2004; Shaya & Olling 2010; Tokovinin & Lépine 2012; Alonso-Floriano et al. 2015, e.g.,). Recent work (Oh et al. 2017; Andrews et al. 2017a; Andrews et al. 2017b; Oelkers et al. 2017; Price-Whelan et al. 2017; El-Badry & Rix 2018; Simpson et al. 2018; Bochanski et al. 2018; El-Badry & Rix 2018; Jiménez-Esteban et al. 2019; Andrews et al. 2019, e.g.,) has shown the power of using Gaia to study co-moving pairs. The availability of 6D phase space information for millions of stars from Gaia DR2 and complementary data from ground-based spectroscopic surveys provide a unique opportunity to understand the nature of co-moving pairs. However, ab initio simulations of galaxy formation currently do not offer the resolution necessary to interpret the small-scale phase space structure observed in Gaia data.

In Kamdar et al. 2019 (hereafter K19) we presented simulations that resolve the dynamical evolution of individual stars comprising the disk over the past 5 Gyr. A key prediction of the fiducial simulation presented in K19 is the high fraction of pairs at large separations (up to $20$ pc) and at low relative velocities (up to $1.5$ km s-1) that were born together (co-natal). In this Letter we use Gaia DR2 and LAMOST DR4 to test the predictions presented in K19 and explore the nature of co-moving pairs.

II. Simulations & Data

In K19 we introduced a new set of simulations that were the first of their kind to model the full population of stars (younger than 5 Gyr) that comprise a Milky Way-like disk galaxy. All stars are born in clusters with a range of initial conditions informed by observations and detailed simulations. The dynamical evolution of 4 billion stars was performed with orbit integration of test particles coupled to a realistic time-varying galactic potential, which includes a disk, halo, bulge, bar, spiral arms, and GMCs (as perturbers). These simulations predict a rich structure in the combined phase and chemical space that should inform our understanding of the nature of clustered star formation.

The fiducial simulation presented in K19 has both non-axisymmetries in the potential (bar, spiral arms, and giant molecular clouds) and clustered star formation. We also ran two control simulations: 1) A simulation with a static axisymmetric potential and with clustered star formation. 2) A simulation with non-axisymmetric perturbations (with bar and spiral arms) but with no clustered star formation (NCSF simulation hereafter). The three different simulations allow us to study the scales at which structure due to clustered star formation and structure due to resonances by the non-axisymmetries in the Milky Way will manifest itself on the combined chemodynamical space.

An accurate error model is essential for comparisons between simulations and observations. The error in the parallax is approximated by a simple error floor of $0.04$ mas (Lindegren et al. 2018). The dependence of proper motion and radial velocity errors is a complex function of several parameters; we fit a gaussian mixture model (GMM) with 20 components to the combined ($G,G_{\rm{BP}}-G_{\rm{RP}},\sigma_{\rm{\mu_{\alpha^{*}}}},\sigma_{\rm{\mu_{\delta}}}$) and ($G,G_{\rm{BP}}-G_{\rm{RP}},\sigma_{\rm{RV}}$) spaces respectively and sample from the conditional distributions for the respective errors given $G$ and $G_{\rm{BP}}-G_{\rm{RP}}$.

We start with the 6D Gaia DR2 (Gaia Collaboration et al. 2018) catalog from Marchetti et al. 2018. We cross-matched with LAMOST DR4 catalog (with duplicates removed) for the same 1 kpc sphere around the solar position as in the simulations. LAMOST DR4 (Deng et al. 2012) provides $T_{\rm eff}$, [Fe/H], radial velocities, and $\log g$ using LAMOST’s spectroscopic pipeline (Wu et al. 2011; Luo et al. 2015). Only stars that have radial velocity measurements in Gaia DR2 are considered for the work presented here since LAMOST RVs have errors between $5-7$ km s-1.

We impose the following selection criteria on the Gaia data considered in this analysis: (1) number of visibility periods $\geq 8$, (2) number of RV transits $\geq 3$, (3) re-normalized unit weight error $\leq 1.6$ (the results presented in this work do not change when making the more conservative $\leq 1.4$ cut used in Lindegren 2018), and (4) bad RVs found in (Boubert et al. 2019) have been removed. We also require SNR${}_{i}>40$ in the LAMOST data to ensure small measurement uncertainties. After making these quality cuts, the quoted mean and median uncertainty on [Fe/H] is 0.037 and 0.024 dex respectively – we emphasize that it is the relative difference between metallicities that is important for this work, rather than the absolute metallicity scale. We have also identified a significant correlation between the metallicity difference and $T_{\rm eff}$ difference for pairs of stars, which we interpret as a systematic uncertainty on the derived metallicities. To limit this systematic uncertainty, we restrict our analysis to pairs with a temperature difference of $\Delta T_{\rm eff}<200$ K.

The overall selection function for Gaia is not critical for this work because we only consider fractional quantities when comparing between data and simulations. We have also explored the impact of the LAMOST selection function on the population of pairs. We find for example no statistically significant difference in the distribution of physical separations for pairs with and without LAMOST data, leading us to conclude that the LAMOST selection function does not result in a biased population of pairs.

Since the simulations in K19 do not include binarity, it is important to discuss how to isolate the phase space signature of disrupting star clusters and avoid contamination from wide binaries in the data. For our analysis, we will only consider pairs with separations greater than $2$ pc, since the Jacobi radius is $\approx 2$ pc for $\sim 1M_{\odot}$ stars in the solar neighborhood (Yoo et al. 2004; Jiang & Tremaine 2010). In addition, Jiang & Tremaine 2010 predict a trough in the distribution of wide binaries between $\sim 1-10$ times the Jacobi radii, which corresponds to separations between $\sim 2-18$ pc in the solar neighborhood. Consequently, we expect the contamination from wide binaries for the spatial scales being considered in this work to be low.

Moreover, co-moving pairs could be a part of known bound clusters, moving groups or OB associations in the Milky Way. We choose to exclude pairs from these to isolate the signal of star clusters that are disrupting. The exploration of the phase space signature in data and simulations for bound objects is left for future work. Similar to Oh et al. 2017, we form an undirected graph where stars are nodes, and edges between the nodes exist for co-moving pairs of stars. Consequently, a star could have multiple co-moving neighbors, and a pair of stars could be directly or indirectly connected via a sequence of edges. The graph is then split into connected components – a connected component is a subgraph of the original graph in which any two nodes are connected to each other by a path – to pick out only mutually exclusive pairs. This additional selection criteria efficiently filters out known open clusters, OB associations, and moving groups.

III. Co-moving Pairs in Simulations & Data

In this section we investigate co-moving pairs in the simulations and the data. We begin with a pair-wise comparison of the fraction of pairs as a function of separation and velocity difference in the simulations and in the data in Figure 1. The left panel of Figure 1 shows the fraction of co-moving pairs as a function of $\Delta r$ for the different simulations (dashed lines of different colors) and the data (black solid line) after imposing a $\Delta v<1.5$ km s-1 cut. At smaller separations, clustering seems to be stronger in the data than in the simulations. The discrepancy at the smallest separations could be due to contamination from wide binaries or hint at complexities in the data regarding star formation that are not adequately modeled in the fiducial simulation.

The right panel of Figure 1 shows the fraction of co-moving pairs as a function of $\Delta v$ for the different simulations and the data with the cut $2<\Delta r<20$ pc. The axisymmetric simulation has a notably higher fraction of pairs at lower $\Delta v$ due to the absence of any large-scale scattering. The NCSF simulation has a notably lower relative fraction at lower $\Delta v$ due to the absence of clustered star formation. The fiducial simulation, which has both clustered star formation and non-axisymmetries in the potential, is in excellent agreement with the data.

Figure 2 shows the co-natal fraction of co-moving pairs in the fiducial simulation presented in K19. The simulation predicts that pairs of stars with a velocity difference up to $1.5$ km s-1 but a high separation of up to $20$ pc are highly likely to be born together. The large physical separation suggests that the co-moving pairs of stars are likely part of a disrupting star cluster – studying such co-natal pairs in data will provide key insights into star formation in the disk. Consequently, using the simulation as a prior of sorts, we will use the selection box $2<\Delta r<20$ pc and $\Delta v<1.5$ km s-1 to both avoid most wide binaries and look for signatures of dissolving star clusters in the data.

While, as we can see in the left panel of Figure 2, the selection box contains a large number of co-moving pairs that are co-natal, there is still some contamination (about $\sim$20%) from field pairs. To test whether we find a similarly high co-natal fraction in data, we can use metallicities since stars born in the same cluster are believed to have essentially identical metallicities (De Silva et al. 2007; Bovy 2016, e.g.,), modulo atomic diffusion (Dotter et al. 2017). The right panel of Figure 2 shows expectations for the metallicity ([Fe/H]) difference distribution of co-moving pairs. The black line shows pairs of stars known to be co-natal in the selection box, the blue line shows all pairs of stars in the selection box, and the red line shows random field pairs. The measurement uncertainty in [Fe/H] in the simulations is $0.03$ dex. The selection criteria are effective at picking out co-natal pairs but there is still enough contamination from field pairs for there to be a tail at larger metallicity differences. The blue line will be used to compare to the Gaia/LAMOST cross-match.

The top panel of Figure 3 shows the metallicity difference ($|\Delta$ [Fe/H]$|$) for co-moving pairs with $2<\Delta r<20$ pc and $\Delta v<1.5$ km s-1. The agreement between the data and the simulation (which has a fiducial $\sigma_{\rm[Fe/H]}=0.03$ dex as a mock observational uncertainty) is very good. Moreover, the distribution of metallicity differences is significantly narrower than the field metallicity difference distribution, suggesting that the co-moving pairs we have identified are co-natal.

The bottom left panel of Figure 3 shows the fraction of pairs with $|\Delta$ [Fe/H]$|<0.1$ as a function of different velocity cuts after already imposing $2<\Delta r<20$ pc. The fraction of stars with a low metallicity difference quickly drops for velocity differences larger than $\sim 1-1.5$ km s-1 in both the data and the simulation. The fraction of pairs with low metallicity differences approaches that of the field pairs’ at progressively larger $\Delta v$ cuts.

The bottom right panel of Figure 3 shows the fraction of pairs with $|\Delta$ [Fe/H]$|<0.1$ as a function of different separation cuts after already imposing $\Delta v<1.5$ km s-1. The fraction of stars with a low metallicity difference drops relatively smoothly as a function of separation for both the simulation and the data – this trend is also clearly visible in Figure 2. The fraction of pairs with low metallicity differences approaches that of the field pairs’ at progressively larger $\Delta r$ cuts. Taken together, the results in Figure 3 indicate that the final selection cuts of $2<\Delta r<20$ pc and $\Delta v<1.5$ km s-1, lead to a fairly clean sample ($\sim 80\%$ pure) of co-moving pairs that were born together.

The simulation offers the opportunity to identify properties of the birth sites of the co-natal pairs. The left panel of Figure 4 shows the mass of the birth cluster for the co-natal co-moving pairs in the selection box in Figure 2 compared to the overall cluster mass function (CMF). The clusters that give birth to co-natal co-moving pairs have notably higher masses relative to the overall CMF. This is due to the much larger number of pairs of stars that are possible with an initially large number of stars in the cluster. The right panel of Figure 4 shows the age of the birth cluster for the co-natal co-moving pairs in the selection box shown in Figure 2 compared to the overall cluster age function. The clusters that give birth to co-natal co-moving pairs are mostly younger than $1$ Gyr. The age distribution of these clusters offers a noteworthy parallel to the discussion of visibility timescale presented in K19, where we argued that the phase space signature of stars born in clusters should be visible up to $1$ Gyr.

IV. Summary

In this Letter we have presented evidence that co-moving pairs of stars (identified as having relative physical separation $2<\Delta r<20$ pc, and relative velocity $\Delta v<1.5$ km s-1) were very likely born together. Moreover, we have shown that the distribution of pair-wise velocities and physical separation is sensitive to both clustered star formation and non-axisymmetries of the Galactic potential. The simulation presented in K19 is able to reproduce both of these distributions. Furthermore, the simulation predicts that the co-natal co-moving pairs were born in preferentially high-mass and relatively young clusters relative to the field population. Table 1 lists all 111 pairs in our final catalog and their properties.

The primary metric used here for determining whether a pair is co-natal is the metallicity difference ($|\Delta$[Fe/H]$|$). For this we used data from the low resolution LAMOST spectroscopic survey since it contains far more stars in the solar neighborhood than other surveys. Follow-up high resolution spectra of these co-moving pairs is required to confirm the chemical similarity of the pairs, which in turn would bolster the argument that the pairs share a common birth site. Gaia DR4 will deliver SNR $\sim 50$ spectra for stars with $G\sim 12$, allowing calculations of [Fe/H] with uncertainties $\leq 0.05$ dex (Recio-Blanco et al. 2016; Ting et al. 2017). The simulations predict that we could study thousands of co-moving pairs due to disrupting star clusters in the solar neighborhood using Gaia alone.

Co-moving pairs offer novel constraints on the nature of clustered star formation and the recent dynamical history of the disk. However, pairs offer a limited ($N=2$) view of the phase space structure of stars in the disk. Additional insight will be gained by considering the general clustering properties of stars in various phase space projections. This will be the subject of future work.

References