LACEwING: A New Moving Group Analysis Code

Within the LACEwING framework, each NYMG is represented by two triaxial ellipsoids described by base values UVW and XYZ, semimajor and semiminor axes ABC and DEF, and Tait-Bryant angles UV, UW, VW, XY, XZ, and YZ, which are used to rotate the ellipses around the W, V, U and Z, Y, X axes, respectively. This is unlike the BANYAN II (Gagné et al. 2014b) rotation order U, V, W (X, Y, Z).

For each NYMG, we use 100,000 Monte Carlo iterations within the triaxial ellipsoids to generate a realistic spread of UVW velocities. Using inverted matrices from Johnson & Soderblom (1987), we convert the UVW velocities into μ_α, μ_δ, and γ for a simulated star at the α and δ of the target, at a standard distance, which we have chosen to be 10 pc in analogy to the absolute magnitude. The lengths of the proper motion vectors are later used to estimate the kinematic distance, which can be compared to a parallax measurement if it exists.

As an example, we use the bona fide β Pic member AO Men (K4 Ve). At the α and δ of AO Men and a distance of 10 pc, an object sharing the mean UVW velocity of β Pic should have μ_α = −45 mas yr⁻¹, μ_δ = 289 mas yr⁻¹, and γ = +16 km s⁻¹.

With the predicted values of μ_α, μ_δ, and γ, and at least some measured kinematic data on the star, the code derives (up to) four of the following metrics:

• 1.  
  Proper motion. The code splits the measured proper motion into components parallel (μ_∥) and perpendicular (μ_⊥) to the predicted proper motion vector; the best match is where the perpendicular component is zero. Our goodness-of-fit metric for proper motions is

  $$\begin{eqnarray}&&{\psi }_{\mu }^{2}=\displaystyle \frac{{\mu }_{\perp \mathrm{star}}^{2}}{{\sigma }_{{\mu }_{\perp \mathrm{star}}}^{2}+{\sigma }_{{\mu }_{\perp \mathrm{pred}}}^{2}}.\end{eqnarray} \tag{ 1 }$$

  If ${\mu }_{\parallel \mathrm{star}}$ has the opposite sign from ${\mu }_{\parallel \mathrm{pred}}$ and is more than 1σ from zero, we add 1000 to ψ_μ to ensure a poor goodness-of-fit metric.For our example object AO Men, the real proper motion (μ_α = −8.3 mas yr⁻¹ and μ_δ = 72.0 mas yr⁻¹) is along a vector very similar to the predicted object, leading to a perpendicular component of 2.7 mas yr⁻¹ and a ψ_μ of 0.05σ.
• 2.   
  Distance. Kinematic distance is derived from a ratio of the star's measured proper motion and the predicted proper motion calculated for a member at 10 pc:

  $$\begin{eqnarray}&&\displaystyle \frac{{\mathrm{Dist}}_{\mathrm{pred}}}{10}=\displaystyle \frac{{\mu }_{\mathrm{para},\mathrm{pred}}}{{\mu }_{\mathrm{para},\mathrm{star}}}.\end{eqnarray} \tag{ 2 }$$

  If a trigonometric parallax distance exists, the resulting goodness-of-fit metric is

  $$\begin{eqnarray}&&{\psi }_{\mathrm{Dist}}^{2}=\displaystyle \frac{{\mathrm{Dist}}_{\mathrm{star}}^{2}-{\mathrm{Dist}}_{\mathrm{pred}}^{2}}{{\sigma }_{\mathrm{Dist},\mathrm{star}}^{2}+{\sigma }_{\mathrm{Dist},\mathrm{pred}}^{2}}.\end{eqnarray} \tag{ 3 }$$

  For our example object AO Men, the two parallel components of the proper motion are 72.4 mas yr⁻¹ (real) and 293.5 mas yr⁻¹ (predicted), which yields a ratio of 4.05 and a distance of 40.5 pc. The actual distance (measured by Hipparcos; van Leeuwen 2007—we are not using the Gaia DR1 results from Gaia Collaboration et al. 2016) is 38.6 pc.
• 3.  
  Radial velocity. The code compares the measured radial velocity to the predicted γ. In this case, the goodness-of-fit metric is

  $$\begin{eqnarray}&&{\psi }_{\gamma }^{2}=\displaystyle \frac{{\gamma }_{\mathrm{star}}^{2}-{\gamma }_{\mathrm{pred}}^{2}}{{\sigma }_{\gamma ,\mathrm{star}}^{2}+{\sigma }_{\gamma ,\mathrm{pred}}^{2}}.\end{eqnarray} \tag{ 4 }$$

  For AO Men, the predicted γ of an ideal member of β Pic is +16.04, while the measured γ of AO Men is +16.25. The difference between the two is 0.1σ.
• 4.  
  Position. With either trigonometric parallax or kinematic distance, the code uses the α, δ and distance to determine how near the star is to the moving group or cluster. As the moving groups are defined with freely oriented ellipses, this requires a matrix rotation to bring a 100,000 element Monte Carlo approximation of the stellar position uncertainty into the coordinate system of the moving group. The goodness-of-fit metric is

  $$\begin{eqnarray}&&{\psi }_{{Pos}}^{2}=\displaystyle \frac{{X}_{\mathrm{star}}^{2}}{{\sigma }_{X}^{2}+{\sigma }_{D}^{2}}+\displaystyle \frac{{Y}_{\mathrm{star}}^{2}}{{\sigma }_{Y}^{2}+{\sigma }_{E}^{2}}+\displaystyle \frac{{Z}_{\mathrm{star}}^{2}}{{\sigma }_{Z}^{2}+{\sigma }_{F}^{2}}.\end{eqnarray} \tag{ 5 }$$

  AO Men and its measured parallax put it 30 pc (2.1σ) from the center of the β Pic moving group.

While the last test for spatial position is not a standard convergence test, it is useful for preventing the code from identifying members of spatially concentrated groups (e.g., Hyades, TW Hya) in physically unreasonable locations.

Each of the four goodness-of-fit metrics has a different characteristic median value when analyzing bona fide members, and that value varies slightly between groups. The median ψ_μ value is particularly small (∼0.02) when matching bona fide members to their correct groups, while ψ_Dist and ψ_γ are much larger (∼0.5), and ψ_Pos was generally around 2. All goodness-of-fit values tended to be significantly larger when matching stars to the wrong groups, although ψ_μ was again the most sensitive discriminator.

At this stage, the LACEwING code has produced up to four different goodness-of-fit metrics estimating the quality of match between a star and one of the NYMGs. We combine these metrics in quadrature:

$$\begin{eqnarray}&&{\psi }^{2}\times {N}_{\mathrm{metrics}}={\psi }_{\mu }^{2}+{\psi }_{\mathrm{Dist}}^{2}+{\psi }_{\gamma }^{2}+{\psi }_{\mathrm{Pos}}^{2}.\end{eqnarray} \tag{ 6 }$$

All tests indicated that the code recovers more members when each metric is given equal weight, as shown here.

We now want to derive functions to transform the goodness-of-fit values into membership probabilities. This is not simple, because the groups overlap with each other to different degrees. Some, like Ursa Major, have distinct UVW velocities that do not overlap with any other group. Others, like the Hyades, have unique and compact XYZ spatial distributions. Most, however, overlap in UVW and XYZ with other moving groups.

The additional complication is that there are seven different possible combinations of data, each with its own goodness-of-fit value that we need to calibrate separately, for each moving group considered:

• 1.  
  α, δ, μ_α and μ_δ.
• 2.  
  α, δ, Dist.
• 3.  
  α, δ, γ.
• 4.  
  α, δ, μ_α and μ_δ, Dist.
• 5.  
  α, δ, μ_α and μ_δ, γ.
• 6.  
  α, δ, Dist, γ.
• 7.  
  α, δ, μ_α and μ_δ, Dist, γ.

To quantify these probabilities, we generate a simulation of the Solar Neighborhood that we run through the LACEwING algorithm. The results of the simulation will be used to determine the relationship between the goodness of fit and membership probability.

When generating the simulation, a random number generator first assigns the new star to one of the groups, in proportion to the population of each group as derived in Section 5.1. Thus, more field stars will be generated than TW Hya members because the field is more populous.

To generate UVW and XYZ values for a simulated star, we use the freely oriented ellipse parameters of the moving groups (from Section 5.1). If the assigned group is one of the NYMGs, a UVW velocity and XYZ position is generated under the assumption that the velocity and spatial distributions are both Gaussian. If the assigned group is the field, the UVW velocities are assumed to be Gaussian distributions, but the XY positions are drawn from a uniform distribution truncated at semimajor and semiminor axes D and E, and the Z position above the plane is drawn from an exponential distribution with a scale height of 300 parsecs (Bochanski et al. 2010) truncated at semiminor axes F.

Once a star has been generated, the UVW and XYZ values are converted into equatorial coordinates α, δ, π, μ_α, μ_δ, and γ. Measurement errors are generated in the form of Gaussian distributions scaled to position to 005, μ to 5 mas yr⁻¹, π to 0.5 mas, and γ to 0.5 km s⁻¹, and added to the equatorial values. The simulated stars are then run through the LACEwING algorithm, comparing them to every moving group using all seven combinations of data. The goodness-of-fit values are then recorded.

To produce the actual relations that yield membership percentages, we then take all of the stars compared to a particular moving group X, using a combination of data Y. The goodness-of-fit values are binned into 0.1 goodness-of-fit value bins. Within each bin, we calculate the fraction of stars in that bin that were actually generated as members ("true" members) of moving group X.

These functions are best described by Gaussian Cumulative Distribution Functions (CDFs) with normalized amplitude. The fraction of stars that were "true members" as a function of goodness-of-fit value (for β Pic) is shown in Figure 1 as histograms, along with the fitted functions. The coefficients of the Gaussian CDFs are saved and used to derive membership probabilities given a goodness-of-fit value.

Zoom In Zoom Out Reset image size

In order to handle the case where a star is already known to be young (and the probability of membership in an NYMG is higher), there is a second calibration of LACEwING. In this mode, a subset containing all the NYMG members and an equally sized portion of the field stars (to account for the young field; see Section 8) is retained, for a 1:1 ratio of NYMG members:young field is selected. This subset is binned into 0.1 goodness-of-fit bins, fractions of stars are calculated in the same way as the above "field star" sample, and a different set of Gaussian CDF coefficients is produced (the "young star" curves for β Pic are shown in Figure 2).

Zoom In Zoom Out Reset image size

These curves are different for different moving groups; the results of Coma Ber in Field Star mode are shown in Figure 3. With 16 groups, 7 combinations of data, and 2 modes of operation, there are 224 sets of coefficients that make up this implementation of LACEwING.

Zoom In Zoom Out Reset image size

For AO Men, our combined goodness-of-fit value was 0.52. In the field star mode histogram (Figure 1) and the case of having complete kinematic data, 3632 of the 8 million stars scored between 0.5 and 0.6, and only 1032 of those (28%) were generated as members of β Pic. If we use the young star mode histograms (Figure 2), 61% of the stars in the 0.5–0.6 bin were genuine members of β Pic. Using the actual curve fits to compute the probabilities of β Pic membership for AO Men yields 26% (field star mode) and 57% (young star mode), with an additional 3% (field star)/13% (young star) chance that it is actually a member of Columba according to those curves.

The membership probabilities given by LACEwING are ultimately the complement of the contamination probability: the probability that a star is a member of a given group and not something else. This is a subtly different question than, "what group X is star Z a member of?" The latter question involves comparing different NYMGs directly, and is much more difficult to answer. When interpreting LACEwING probabilities, it is important to keep in mind that LACEwING does not force all membership probabilities to add up to 100% in the way that BANYAN and BANYAN II do; each of the probabilities of membership is an independent assessment of "Group X or not Group X?" connected only by the fact that all the probability coefficients are derived from the same simulation. Probabilities indicated by LACEwING may add up to more than 100% if the uncertainties on the input parameters are larger than the typical values in the simulation. In practice, taking the NYMG that is matched with the highest membership probability is an excellent means of identifying memberships.

It is important to generate enough simulated stars to properly sample the membership function for even the smallest group—a minimum of roughly 1000 stars are necessary to fit the Gaussian CDF correctly. For the particular implementation of LACEwING presented in this paper, the smallest group was η Cha (6 members, out of 40,902), which required at least 6.8 million simulated stars; 8 million were calculated.

From Figure 1, in field star mode and with all kinematic data, even the best possible match to β Pic (goodness-of-fit = 0) has only a 92% probability of membership in β Pic. That is, 8% of objects that match β Pic perfectly are not members. With only proper motion, the maximum possible probability of membership in β Pic is only 5%. The young star mode yields higher probabilities; Figure 2 shows a maximum probability of 100% for β Pic members with all information, and a maximum membership probability of 23% for only proper motion. If we plot the cumulative number of β Pic members recovered as a function of minimum probability accepted (Figure 4), we see that we have to set a minimum threshold of 10% membership probability to recover 90% of the members in the best possible case where all kinematic data are available.

Zoom In Zoom Out Reset image size

Setting a low membership probability cutoff means a much larger false positive contamination, as shown in Figure 5. Selecting an adequate cutoff requires balancing the recovery rate with the contamination rate (Figure 6; similar figures describing the BANYAN II Bayesian models can be found in Figures 5 and 6 of Gagné et al. 2014b). The false positive rates highlight the danger of using kinematics alone to identify young stars: any kinematically selected survey needs to use other spectroscopic and photometric youth indicators to weed out false positives.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Rough guidelines used throughout the rest of this work, and in the stand-alone version of LACEwING, are that probabilities of 66% and higher are high probabilities of membership, 40%–66% and above are moderate probability, 20%–40% are low probability, and below 20% is too low to consider meaningful.

Technical details on using LACEwING, recalibrating LACEwING, and incorporating it into other codes are given in Appendix A.

The TRACEwING code uses an epicyclic approximation of Galactic orbital motion (Makarov et al. 2004) to trace the positions of stars back in time using their current measured motion, in increments of 0.1 Myr. It compares the positions of single objects to an NYMG, which is represented by stored freely oriented ellipse parameters fit using the same process used in Section 2.1. Based on the equal-volume radius of the group, the TRACEwING code presents a single-valued representation of how close the target is to the moving group as a function of time.

TRACEwING is essentially two separate steps, carried out by two different programs:

• 1.  
  A program that uses the epicyclic kinematic approximation to trace all bona fide members of an NYMG back in time, and fits freely oriented ellipsoids to the ensemble at each time step, saving the parameters for future use.
• 2.  
  A program that uses the epicyclic kinematic approximation to trace a single star back in time, and compares its positions to saved moving group ellipsoids at each time step.

With epicyclic traceback, the effects of Galactic orbital motion are approximated by use of sine and cosine functions, controlled by Oort constants and a vertical oscillation parameter. For TRACEwING, we use the equations of position (relative to the Sun, as a function of time in Myr) given in Makarov et al. (2004) and reproduced here:

$$\begin{eqnarray}&&X={X}_{0}+{U}_{0}{\kappa }^{-1}\sin \kappa T+({V}_{0}-2{{AX}}_{0})(1-\cos \kappa T){(2B)}^{-1},\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}Y & = & {Y}_{0}-{U}_{0}(1-\cos \kappa T){(2B)}^{-1}\\ & & +{V}_{0}(A\kappa T-(A-B)\sin \kappa T){(\kappa B)}^{-1}\\ & & -2{X}_{0}A(A-B)(\kappa T-\sin \kappa T){(\kappa B)}^{-1},\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&Z={Z}_{0}\cos \nu T+{W}_{0}{\nu }^{-1}\sin \nu T.\end{eqnarray} \tag{ 9 }$$

In the above equations, A and B are the Oort constants (from Bobylev 2010; A = +0.0178 km s⁻¹ pc⁻¹, B = −0.0132 km s⁻¹ pc⁻¹); κ is the planar oscillation frequency, $\sqrt{-4B\times (A-B)}$; and ν is the vertical oscillation frequency, 2π/85 Myr⁻¹ (Makarov et al. 2004). The approximation deviates noticeably from an unperturbed linear traceback motion after 10 million years.

For the first step of TRACEwING, we calculate 4000 Monte Carlo tracebacks for each bona fide member, in time increments of 0.1 Myr. At each time step, freely oriented ellipses are fit to the positions of the members of the moving group for each of the 4000 Monte Carlo trials, and averaged to produce a mean position, extent, and orientation of the group at that time. These parameters are saved for comparison to individual stars.

For the second step, we take a target object of interest, which may be any object with full kinematic information—star, brown dwarf, or planet. To determine the potential memberships of the target object, we generate and trace back 20,000 trials within each of the 1σ, 2σ, and 3σ uncertainties on their observational (equatorial) positions and motions. At each time step, the distance between the object and the previously calculated position of the moving group (both in Cartesian Galactic XYZ coordinates) is calculated, and an equal-volume radius ($\tfrac{4}{3}\pi {r}^{3}=\tfrac{4}{3}\pi {abc}$) is used as the effective radius of the group. The targets are then visually classified by whether the 1σ, 2σ, or 3σ positions potentially place them within the effective radius of the group at the time of formation. The traceback of the bona fide β Pic member AO Men is displayed in Figure 7, and shows a star that was plausibly within the confines of β Pic at the time of formation.

Zoom In Zoom Out Reset image size

Epicyclic approximations do not take into account the gravitational influence of other stars, molecular clouds, Gould's Belt, or the Galactic disk and bar itself. This is most pronounced in open clusters, where the stars themselves are gravitationally bound to each other, but sets limits on the reliability of the technique for moving groups as well.

To quantify the limits of the technique, we perform two tests. First, we simulate a "real" moving group of stars, move it forward in time, "observe" it, and trace it back in time to see what a genuine NYMG of various ages should look like traced back to its origin. Second, we select unrelated field stars with a velocity and spatial distribution similar to known NYMGs, and move it back in time as a "fake" NYMG. The upper limit on reliability of the traceback technique is the point at which the "real" moving group is indistinguishable in volume from the "fake" one.

Our "real" moving group is a simulated group of 50 stars with a Gaussian velocity dispersion of 1.5 km s⁻¹ (Preibisch & Mamajek 2008) typical of the Sco-Cen star-forming region, and a uniform spatial distribution with a radius of 5 pc. These values are perhaps smaller than most clusters, but provide a best-case scenario. The 50 stars were moved forward in time in steps of 0.1 Myr using the epicyclic approximation and their positions (with the mean position subtracted off, so all stars were near the origin at the final epoch) were saved at 5, 8, 12, 25, 45, 50, 125, 250, 400, and 500 Myr intervals. To observe the stars, we generated new UVW velocities using the individual stars' change of position over the last 0.1 Myr of the simulated time range, e.g.,

$$\begin{eqnarray}&&U=0.9778\displaystyle \frac{\mathrm{km}\,\mathrm{Myr}}{{{\rm{s}}}^{-1}\,\mathrm{pc}}{\rm{\Delta }}X\displaystyle \frac{\mathrm{pc}}{0.1\,\mathrm{Myr}}.\end{eqnarray} \tag{ 10 }$$

We then converted all UVWXYZ values to equatorial coordinates, and applied randomly generated "observational" errors: 0.5 mas π uncertainties, 10 mas yr⁻¹ μ uncertainties in each axis, and 1 km s⁻¹ γ uncertainties. These collections of stars were run back in time the same way as before to determine the apparent size at formation. To test the trivial case (perfect information), we added no observational errors to the generated cluster of stars and traced them back in time.

To assemble the group of field stars, we searched the Extended Hipparcos catalog of Anderson & Francis (2012) to find stars with μ, π, and γ distributed according to the present-day median velocity dispersion and spatial distribution parameters of our unbound moving groups, with a velocity dispersion of 1.6 km s⁻¹ and spatial distribution within 11.5 pc. Our fake moving group is centered on XYZ = (−5, −5, 20) pc and UVW = (−5, −5, −5) km s⁻¹, a well-populated region of velocity space not populated by any known NYMG. The 15 selected stars are given in Table 2.

In the trivial case with no observational uncertainties (Figure 8), the synthetic group traces back to having its smallest volume (although it is larger than the 5 pc initial radius) near the actual time of formation. In the more realistic case (Figure 9), the large measurement errors mean that the minimum size of the simulated moving group appears to be right now, and the apparent effective radius of the simulated moving group at time of formation is 40 pc. In Figure 10 we plot the estimated volume-at-formation of the simulated moving group at various ages, and the fake moving group of field stars moved back in time.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the case of the fake moving group made up of unrelated field stars (Figure 11), the volume is smaller than the simulated moving group (Figure 10) after 125 Myr. The epicyclic approximation's Galactic shear causes the simulated group to have a larger present-day velocity dispersion than is currently expected for the known NYMGs, suggesting that we are missing outlying members of the known groups. Crucially, these stars should trace back toward the origin of the moving group as they do in our simulated moving group example, unlike the outliers we intend to remove from current membership lists.

Zoom In Zoom Out Reset image size

The Gaia mission will provide astrometry that is several magnitudes more precise than is currently available for stars brighter than G ≈ 20.7. To determine the effectiveness of Gaia data, we have repeated our first test using measurement uncertainties expected for G ∼ 15 sources (Gaia Collaboration et al. 2016) of 50 μas for parallaxes, 40 μas for positions, and 25 μas yr⁻¹ for each axis of proper motions. The simulated Gaia-precision data demonstrate significant improvements. The size of the simulated burst of stars shrinks by a factor of 1000 (roughly linearly with the increased precision). As late as 25 Myr (Figure 12), there is still a minimum (although not the correct radius) in the spatial distribution of the stars near the actual age of the group, and the radius is now relatively constant over time.

Zoom In Zoom Out Reset image size

We can conclude a few things about the efficacy of tracebacks from these tests, presented here in no particular order:

• 1.  
  The effectiveness of tracebacks is almost entirely limited by the current measurement precision of the parameters α, δ, π, μ, and γ. After roughly 125 Myr, the cumulative effects of measurement uncertainties make it impossible to distinguish between a group of stars spreading out from a single point of origin and an unphysical collection of field stars on roughly parallel tracks. We cannot, therefore, comment on the existence of any moving groups older than 125 Myr based on epicyclic tracebacks alone. As measurement uncertainties shrink, we will gradually approach the "perfect" case (Figure 8) where genuinely related stars will trace back to volumes smaller than unrelated stars selected with a similar velocity dispersion.
• 2.   
  We cannot say anything about the membership of objects that trace back to within the boundary of the moving group at the time of formation—with the currently available measurement precision they could be true members or coincidental field stars. However, stars that do not trace back to possibly be within the confines are a different case. We found that only 6% of simulated members in our 45 Myr old sample were not plausibly within the boundary of the moving group (at the 1σ positional uncertainty level) 45 Myr in the past. In contrast, far more than 6% of actual moving group members (Section 5.1) did not trace back to within the confines of their purported moving group, suggesting that the objects are not real members of the group. These are easily identifiable nonmembers, and as data precision improves, we expect to find more of these objects.
• 3.  
  Very few of the stars in the fake moving group are consistent with possibly being in the center of the moving group (or really, anywhere except the edges of the group), while in the simulated moving group, nearly all the members are consistent with being in the center of the NYMG. This would seem to be one difference between an actual NYMG and an unphysical selection of stars.

The catalog is supplemented by additional data sources drawn from large surveys. This aids in providing useful information about the stars and strengthens the consistency of the source data.

4.1.1. Astrometry

Positions and proper motions were preferentially sourced from the following ICRS catalogs tied to the Hipparcos reference frame:

• 1.  
  Hipparcos 2 (van Leeuwen 2007), 1915 objects, 0.1–3.0 mas position precision, 0.1–5.0 mas yr⁻¹ μ precision.
• 2.  
  UCAC4 (Zacharias et al. 2013), 1527 objects, 10–100 mas αδ precision, 1–10 mas yr⁻¹ μ precision.
• 3.  
  Tycho-2 (Høg et al. 2000), 906 objects, 10–100 mas αδ precision, 1–10 mas yr⁻¹ μ precision.
• 4.  
  PPMXL (Röser et al. 2010), 712 objects, 50–100 mas αδ precision, 5–20 mas yr⁻¹ μ precision.
• 5.  
  2MASS and 2MASS-6X (Skrutskie et al. 2006). 286 objects, 60–120 mas αδ precision, no μ.
• 6.  
  SDSS DR9 (Ahn et al. 2012) 4 objects, 80–200 mas αδ precision, no μ.
• 7.  
  Source papers (from Table 4). These mostly provide proper motions (5–20 mas yr⁻¹ μ precision) for 261 stars found in 2MASS or 2MASS-6X. The other 23 have proper motions copied from a primary star.

There are still 43 objects (35 of which are in the Pleiades) that have no reported proper motions. Parallaxes were sourced from papers listed in Table 5. Where available, parallaxes for objects in multiple systems and multiple observations of the same stars were combined into weighted mean system parallaxes, under the assumptions that the published uncertainties are accurate and that the parallax of every component of the system is the same to within measurement errors.

Table 4.  Source Papers

Citation	Groups Includeda
Eggen (1991)	IC 2391 Supercluster
Barrado y Navascues (1998)	Castor
Makarov & Urban (2000)	Car-Vela
van den Ancker et al. (2000)	Capricornus
Montes et al. (2001b)	Multiple
Montes et al. (2001a)	Multiple
King et al. (2003)	Ursa Major
Ribas (2003)	Castor
Zuckerman & Song (2004)	Multiple
Mamajek (2005)	TW Hya
Casewell et al. (2006)	Coma Ber
López-Santiago et al. (2006)	Multiple
Moór et al. (2006)	Multiple
Zuckerman et al. (2006)	Car-Near
Guenther et al. (2007)	Multiple
Kraus & Hillenbrand (2007)	Coma Ber
Makarov (2007)	Multiple
Platais et al. (2007)	IC 2391 Cluster
Stauffer et al. (2007)	Pleiades
Kirkpatrick et al. (2008)	TW Hya
Mentuch et al. (2008)	Multiple
Teixeira et al. (2008)	TW Hya
Torres et al. (2008, p. 757)	Multiple
da Silva et al. (2009)	Multiple
Guillout et al. (2009)	Other
Lépine & Simon (2009)	β Pic
Shkolnik et al. (2009)	Other
Teixeira et al. (2009)	β Pic, TW Hya
Caballero (2010)	Castor
López-Santiago et al. (2010b)	Multiple
López-Santiago et al. (2010a)	η Cha
Maldonado et al. (2010)	Multiple
Murphy et al. (2010)	η Cha
Nakajima et al. (2010)	Multiple
Rice et al. (2010)	β Pic
Schlieder et al. (2010)	β Pic, AB Dor
Yee & Jensen (2010)	β Pic
Kiss et al. (2011)	Multiple
Riedel et al. (2011)	Argus
Rodriguez et al. (2011)	TW Hya, Sco-Cen
Shkolnik et al. (2011)	TW Hya
Wahhaj et al. (2011)	AB Dor
Zuckerman et al. (2011)	Multiple
Kastner et al. (2012)	Cha
McCarthy & White (2012)	β Pic, AB Dor
Mishenina et al. (2012)	Other
Murphy (2012)	Other
Nakajima & Morino (2012)	Multiple
Schlieder et al. (2012b)	β Pic, AB Dor
Schlieder et al. (2012c)	β Pic, AB Dor
Shkolnik et al. (2012)	Multiple
Schneider et al. (2012a)	TW Hya
Schneider et al. (2012b)	TW Hya
Xing & Xing (2012)	Other
Barenfeld et al. (2013)	AB Dor
Delorme et al. (2013)	Tuc-Hor
De Silva et al. (2013)	Argus
Eisenbeiss et al. (2013)	Her-Lyr
Hinkley et al. (2013)	Other
Liu et al. (2013b)	β Pic
Malo et al. (2013)	Multiple
Mamajek et al. (2013)	Other
Moór et al. (2013)	Multiple
Murphy et al. (2013)	η Cha, Cha
Rodriguez et al. (2013)	Multiple
Schneider et al. (2013)	Other
Zuckerman et al. (2013)	Oct-Near
Weinberger et al. (2013)	TW Hya
Binks & Jeffries (2014)	β Pic
Casewell et al. (2014)	Coma Ber, Hyades
Klutsch et al. (2014)	Multiple
Kraus et al. (2014a)	Tuc-Hor
Gagné et al. (2014b)	Multiple
Gagné et al. (2014a)	TW Hya
Gagné et al. (2014c)	Argus
Malo et al. (2014a)	Multiple
Malo et al. (2014b)	Multiple
Mamajek & Bell (2014)	β Pic
McCarthy & Wilhelm (2014)	AB Dor
Riedel et al. (2014)	Multiple
Rodriguez et al. (2014)	β Pic
Schneider et al. (2014)	Other
Gagné et al. (2015)	Multiple
Murphy & Lawson (2015)	Octans
E. E. Mamajek (2016, private communication)	Multiple

Notes. The papers that make up the membership of the catalog of young stars, along with the groups considered.

^aPapers marked "Multiple" consider multiple groups; papers marked "Other" consider nearby young stars but do not identify them as members of any groups.

Download table as:  ASCIITypeset images: 1 2

4.1.2. Radial Velocity

Radial velocities were assumed to apply to all companions within roughly 2''. Given that it is both possible and likely that different members of a multiple star system have different radial velocities, when different components have independently measured γ, they were combined using a weighted standard deviation to produce a systemic velocity.

Attempts have been made to reduce the double-counting of γ, particularly in cases where a later paper cited a γ taken from one of the catalogs included here. This has particularly been a problem for Barbier-Brossat & Figon (2000) and Kharchenko et al. (2007), which both contain the General Catalog of Radial Velocities (Wilson 1953) where uncertainties were reported as letter codes. Barbier-Brossat & Figon (2000) and Kharchenko et al. (2007) recommend a translation of the letter code into km s⁻¹ uncertainty different from the VizieR version of Wilson (1953; and papers that cite it directly, such as Malo et al. 2013), which it does, leading to nearly identical γ appearing in different sources.

In many cases, γ have been published without uncertainties. Because our weighted standard deviations require an uncertainty, we have invented them where necessary, and flagged them with "e" in our source tables. Radial velocities with quality codes, originating from Wilson (1953), had uncertainties assigned according to the quality notes as suggested in the table notes; where no quality code was available, we have set the errors to 3.7 km s⁻¹, equivalent to letter code "C." Most other papers were given 1 km s⁻¹ uncertainties, as per Table 6.

Table 5.  Survey Papers

Citation	Data Used
Houk & Cowley (1975)	Spectral Types
Houk (1978)	Spectral Types
Houk (1982)	Spectral Types
Andersen et al. (1991)	Multiplicity
Gliese & Jahreiß (1991)	Spectral Types
Hoffleit & Jaschek (1991)	Catalog Names
Kirkpatrick et al. (1991)	Spectral Types
Cannon & Pickering (1993)	Catalog Names
Gatewood et al. (1993)	π
Gatewood (1995)	π
Gatewood & de Jonge (1995)	π
Reid et al. (1995)	Spectral Types
van Altena et al. (1995) (YPC4)	π, Opt phot.
Covino et al. (1997)	γ, Li
Benedict et al. (1999)	π
Söderhjelm (1999)	πa
Voges et al. (1999) (RASS-BSC)	X-ray phot.
Weis et al. (1999)	π
Barbier-Brossat & Figon (2000) (GCMRV)	γ
Benedict et al. (2000)	π
Ducati et al. (2001)	Opt phot.
Høg et al. (2000) (TYCHO-2)	pos., μ, Opt phot.
Voges et al. (2000) (RASS-FSC)	X-ray phot.
Dahn et al. (2002)	π
Gizis et al. (2002)	γ
Henry et al. (2002)	Spectral Types
Nidever et al. (2002)	γ
Torres & Ribas (2002)	πa
Cutri et al. (2003) (2MASS)	pos., NIR phot.
Song et al. (2003)	γ
Thorstensen & Kirkpatrick (2003)	π
McArthur et al. (2004)	π
Pourbaix et al. (2004) (SB9)	Multiplicity
Vrba et al. (2004)	π
Costa et al. (2005)	π
Jao et al. (2005)	π
Lépine & Shara (2005)	μ
Soderblom et al. (2005)	π
Valenti & Fischer (2005)	γ, $\mathrm{vsini}$
Benedict et al. (2006)	π
Gontcharov (2006)	γ
Gray et al. (2006)	Spectral Types
Henry et al. (2006)	π, Opt phot.
Torres et al. (2006)	γ, Hα, Li, $\mathrm{vsini}$
Biller & Close (2007)	π
Close et al. (2007)	Spectral Types
Daemgen et al. (2007)	Multiplicity
Gizis et al. (2007)	π, phot.
Kharchenko et al. (2007)	γ
Scholz et al. (2007)	γ
van Leeuwen (2007) (Hipparcos-2)	π, pos., μ
Ducourant et al. (2008)	π
Fernández et al. (2008)	γ
Jameson et al. (2008)	μ
Gatewood & Coban (2009)	π
Subasavage et al. (2009)	π
Teixeira et al. (2009)	π
Bergfors et al. (2010)	Multiplicity
Blake et al. (2010)	γ
Raghavan et al. (2010)	Multiplicity
Riedel et al. (2010)	π, Opt phot.
Röser et al. (2010) (PPMXL)	position, μ, NIR phot.
Shkolnik et al. (2010)	Multiplicity, γ
Smart et al. (2010)	π
Stauffer et al. (2010)	Catalog Names
Bianchi et al. (2011) (GALEX DR5)	UV phot.
Girard et al. (2011) (SPM4)	Opt phot.
Messina et al. (2011)	Spectral Types
Moór et al. (2011)	γ
von Braun et al. (2011)	π
Ahn et al. (2012) (SDSS DR9)	pos., SDSS phot.
Allen et al. (2012)	μ
Bailey et al. (2012)	γ
Bowler et al. (2012a)	Multiplicity
Bowler et al. (2012b)	Multiplicity
Dupuy & Liu (2012)	π
Faherty et al. (2012)	π
Janson et al. (2012)	Multiplicity
Zacharias et al. (2013) (UCAC4)	pos., μ, Opt, SDSS, NIR phot.
Bowler et al. (2013)	Multiplicity
Cutri et al. (2013) (ALLWISE)	MIR phot.
Kordopatis et al. (2013) (RAVE DR4)	γ
Liu et al. (2013a)	π
Marocco et al. (2013)	π
Dieterich et al. (2014)	π, Opt phot.
Dittmann et al. (2014)	π
Ducourant et al. (2014)	μ, π
Lurie et al. (2014)	π, Opt phot.
Naud et al. (2014)	Multiplicity
Zapatero Osorio et al. (2014)	π
Elliott et al. (2015)	γ
Henden et al. (2015) (APASS DR9)	Opt, SDSS phot.
Mason et al. (2015) (WDS)	Multiplicity
Faherty et al. (2016)	γ

Notes. Additional data sources used in the catalog.

^aParallax replaces Hipparcos data.

Download table as:  ASCIITypeset images: 1 2

4.1.3. Photometry

4.1.4. Multiplicity

LACEwING requires kinematic models of the NYMGs. These are UVW and XYZ ellipsoids fit to genuine members of the groups. To create a list of bona fide members, we have pulled previously identified bona fide members from the Catalog of Suspected Nearby Young Stars (Section 4). We then filtered out probable interlopers from the samples using the TRACEwING code (Section 3). The resulting filtered samples of bona fide moving group members were then fit with ellipsoids.

5.1.1. Initial Member Data

5.1.2. Bona Fide Candidate Filtering

The bona fide list of members was filtered using TRACEwING (Section 3) to identify and remove outliers: stars that could not possibly have been in the same location as the rest of the group at the time of formation. Moving groups were generated from the bona fide list (Figure 14), and then every member of the moving group was traced back to the rest of the group individually. Outliers were defined as being more than 2σ from the location of the NYMG at all times between the minimum and maximum reasonable age for the group, as collected in Table 1. As an example, Figure 15 shows the Tuc-Hor bona fide member CPD-64 17 within the confines of Tuc-Hor over the entire range of quoted ages (30–45 Myr), while Figure 16 shows the Tuc-Hor nonmember HIP 104308 nowhere near Tuc-Hor at any time. After the end of a filtering step, the moving group was recalculated with the refined member list. The process of filtering outliers was repeated until the moving group was self-consistent, which took three or fewer iterations. This reduced the bona fide list to 297 systems.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The volumes of the NYMGs near their reported times of formation are shown in Figure 17, overplotted on Figure 10. With currently available data, the final volume of β Pic is actually smaller than calculated for our synthetic 25 Myr old cluster (Section 3.3), which suggests that it is consistent with being a genuine product of a single burst of star formation. The same is true of AB Dor (125 Myr), suggesting that it too is consistent with being a real moving group, although AB Dor is also close to indistinguishable from the fake cluster of field stars (Section 3.3). The traceback results for AB Doradus in McCarthy & Wilhelm (2014) are similar in implied volume to the TRACEwING traceback of AB Doradus in Figure 14, despite using different epicyclic parameters. This again suggests that the limiting factor in both cases is data precision.

Zoom In Zoom Out Reset image size

The complete table of bona fide members (including discarded nonmembers) is given in the Catalog of Suspected Nearby Young Stars (Section 4), where they are flagged in column "Bonafide" with "B" (for bona fide system primaries), "R" (for rejected system primaries), or "X" (for bona fide systems without sufficient information that could not be filtered with tracebacks or used to construct kinematic models); lower case "b," "r," and "x" flags indicate companions to the respective systems.

As outlined in Section 2, the LACEwING code relies upon triaxial ellipsoid representations of the NYMGs, and an assessment of the population size. Most of the moving group ellipses were created by fitting the filtered selection of bona fide stars. The fitting routine assumes that the groups are triaxial ellipsoids with orthogonal axes. The routine finds the UV plane angle with a linear fit to the projected UV data, de-rotates the data to align that axis with the Cartesian plane, and then repeats the process for the UW plane angle and the VW plane angle. Standard deviations are fit to the de-rotated data and are taken as the axis dimensions of the ellipse. This process is repeated for 10,000 Monte Carlo iterations. The final ellipse parameters are the average of this process, and are shown in Table 7. Two-dimensional projections of the ellipsoids are shown in Figure 18.

Zoom In Zoom Out Reset image size

5.2.1. Groups Whose Properties Did Not Come from Ellipse Fitting

5.2.2. The Field Population

As shown in the top panels of Figure 19, the kinematics of the solar neighborhood as plotted from the XHIP catalog (Anderson & Francis 2012) have a complex structure. Of the NYMGs and open clusters, only the Hyades is readily visible; the remainders of the structures are thought to be the result of Galactic resonances.

Zoom In Zoom Out Reset image size

We have replicated the structures with a by-eye fit of seven ellipsoid components, which correspond to (in terms of the Skuljan et al. 1999 groups) Sirius, Coma Berenices (leading), Coma Berenices (trailing), Pleiades (leading), Hyades (trailing), ζ Hercules, and a broad generic field population. Following Skuljan et al. (1999), all groups are inclined by 25° to the coordinate axes, with the exception of the Pleiades branch at −25° and the unrotated field population.

The bottom panels of Figure 19 demonstrate our synthetic field population, where the top panels plot stars with <50% parallax uncertainties and γ measurements with uncertainties from the XHIP catalog. Figure 20 shows the same plot for the XYZ population, which we have modeled in X and Y as a uniform distribution truncated at a distance of 120 pc (to accommodate groups like Octans and Cha that extend beyond 100 pc), and in Z as an exponential with a scale height of 300 pc, again truncated at a distance of 120 pc.

Zoom In Zoom Out Reset image size

In order to provide an appropriate simulation of members, we must also consider the relative populations of the groups. For these purposes we returned to the Catalog of Suspected Nearby Young Stars (Section 4) and considered all members of the groups, beyond just the groups that survived our bona fide vetting process above.

There are two major sources of incompleteness that must be considered here. First, the more recently discovered groups and older groups (where stars are less obviously youthful) have not been searched as completely for low-mass stars as younger and longer-known groups. For example, none of the 14 known members of the rarely studied χ⁰¹ For group (∼500 Myr) have M dwarf primaries, but 31 of the 38 known members of TW Hya (∼10 Myr) have M dwarf primaries. Second, the continued reliance on the highly magnitude-limited Hipparcos and Tycho-2 catalogs makes it hard to identify members of more distant groups.

To account for this bias, we count only the B, A, F, G, and K primaries (and evolved versions of the same) when tallying membership in these groups. This relies upon three further assumptions. First, that all the groups are similarly complete down to spectral type K7; this will cause more rarely studied groups like 32 Ori and χ⁰¹ For to be underpredicted. Second, that the initial mass function of all of these groups has the same slope; this will overpredict membership in groups with known top-heavy mass functions like the Hyades (Goldman et al. 2013), and in theory underpredict membership in any group with a bottom-heavy mass function. Third, that a K7 spectral type refers to the same mass at the age of every group. This is harder to quantify given the differences in evolutionary tracks between different stellar evolution codes, but should cause younger groups (whose members will be cooler) like TW Hya and β Pic to be slightly underpredicted relative to older groups.

For the field star population, we analyze the Catalog of Suspected Nearby Young Stars for lithium-rich (see Section 5.4) and bona fide members, and find that 271 systems that are bona fide or lithium rich are predicted to be within 25 pc of the Sun. Given the density of star systems within 5 pc (52 systems) there should be 6500 systems within 25 pc of the Sun (a sample which is highly incomplete, with only 2184 systems currently known according to T. J. Henry et al. 2016, in preparation), a ratio of just under 24:1. We therefore take the ratio of young star systems to field stars to be 25:1. Given the further result (see Section 7.23) that less than half of the young stars are in moving groups, we take the ratio of young field stars to moving group and open cluster members to be 1:1. The ratio of all field stars (young and old) to moving group members is therefore 50:1. The population of field stars (split across the seven subgroups, Section 5.2.2) is thus set at 50 times the combined number of all moving group and open cluster members.

As a proof of concept, we use LACEwING on a sample of stars with detectable lithium. Given that kinematics and activity alone are not sufficient to say that an object is young, our process of examining only spectroscopically young objects reduces the possibility that we are attempting to reproduce erroneous membership assignments.

The presence of lithium (specifically the red-optical 6708 Å doublet) is one of the most reliable spectroscopic methods of identifying young stars (Soderblom 2010). Lithium is fused at lower temperatures than hydrogen, and nearly all of it is primordial. Any amount of lithium present in a stellar spectrum (unless the object is an asymptotic giant branch star) is a sign that the object has not fused it yet, either because the object is very young, or because it is a <60M_Jup brown dwarf that never reaches sufficiently high temperatures (Rebolo et al. 1992). The exact age at which lithium is depleted varies based on the mass of the object (Yee & Jensen 2010) and spans ages typical of the NYMGs. In Figure 21, we plot all lithium measurements from the Catalog of Suspected Nearby Young Stars, with fifteen-element moving averages tracing out the lithium depletion as a function of age. At the age of Cha, barely any lithium is gone, while by the age of AB Dor, there is a very clear lithium depletion.

Zoom In Zoom Out Reset image size

Although our survey of papers reporting lithium is not complete, 1877 of the 5350 stars in our catalog of suspected young stars (Section 4) have at least one attempt to measure their Li i λ6707.8 Å EW doublet (EW(Li)), as shown in Figure 21. EW(Li) is frequently reported without uncertainties on the measurements. We have created a lithium sample by selecting objects with EW(Li) > max(10 mÅ, 2σ_Li) from the Catalog of Suspected Nearby Young Stars (Section 4). Where no uncertainty is quoted, the uncertainty is set at 50 mÅ. This matches the limits of several major papers reporting lithium, including Eisenbeiss et al. (2013), Kraus et al. (2014b), and Rodriguez et al. (2013); the largest upper limits reported by Malo et al. (2013) and Moór et al. (2013) are 46 mÅ.

Objects that are companions (according to the WDS and other resources in Section 4) to stars in the lithium sample are also included, yielding a total of 1179 stars. From there, we considered only the primaries of the systems that have at least one lithium-detected member (including 15 primaries that do not themselves have measured lithium) so as not to count star systems more than once, resulting in a list of 1037 star systems.

Objects that were members of groups not within 100 pc were removed. This includes the members of IC 2391 identified by Torres et al. (2008, p. 757) as members of Argus, and a wide variety of objects from Guenther et al. (2007) that belong to more distant regions (e.g., Chamæleon, ρ Ophiuchi, Scorpius–Centaurus). This reduced the size of the lithium sample to 930 systems.

This sample substantially overlaps with the bona fide sample (Section 5.1): 152 systems are in both lists. Some genuinely young stars at the fully convective boundary around spectral type M3 V are excluded from this sample because lithium is fully depleted at their ages, while some more massive stars still maintain lithium even at old ages (e.g., α Cen AB; Mishenina et al. 2012). The sample can be found in the "LiSample" column of the Catalog of Suspected Nearby Young Stars, where system primaries that qualify for the lithium sample are flagged with "L," system primaries that qualify for the lithium sample but have membership in a more distant group are flagged with "F," and primaries that do not have lithium but have a companion that qualifies are flagged with "A." As with the bona fide sample, lowercase letters follow the same rules but denote companions.

6.1.1. Public Codes

6.1.2. Other Codes

The convergence method used in Montes et al. (2001b) is an implementation of techniques used by Eggen (1995). It uses the velocity contributions V_tangential = 4.74 μπ⁻¹ and V_radial to determine the total velocity V_Total. Using the angular separation between the star and the convergent point (λ), it splits V_Total into velocity components parallel (V_T) and perpendicular to (PV) the convergent vector. It also splits the parallel velocity vector V_T into an RV component ρ_C. The two metrics of membership reported are whether PV < 0.1V_T and whether ρ_C and V_radial differ by less than 4–8 km s⁻¹—no percentages are calculated. As used in Montes et al. (2001b), the code considers membership in the Local Association, Hyades Supercluster, IC 2391 Supercluster, Ursa Major moving group, and Castor moving group.

The convergent point method used by Mamajek (2005) is much more similar to the Rodriguez convergence code.

The codes used by Lépine & Simon (2009), Schlieder et al. (2010, 2012b, 2012c), and Kraus et al. (2014a) are quite similar to LACEwING. They take the UVW velocity of a moving group and determine the expected proper motion and radial velocity of a member at those equatorial coordinates (at an implied distance of 1 pc); the proper motion vector is then used to determine a kinematic distance and kinematic RV. These kinematic distances are compared to spectroscopic, photometric, or trigonometric distances (Lépine & Simon 2009), or with spectrophotometric distances and a distance modulus cutoff to enforce the stars falling on the approximately correct isochrone (Kraus et al. 2014a). They lack the XYZ position metric used by LACEwING, and do not convert the proper motion matches or distance matches to membership probabilities. Lépine & Simon (2009) only reports members of β Pic; Schlieder et al. (2010, 2012b, 2012c) report β Pic and AB Dor; Kraus et al. (2014a) only considers Tuc-Hor. It is not clear if the codes consider more than one or two groups in their studies.

The moving group identification techniques used by Shkolnik et al. (2012), Riedel et al. (2014; see also Baines et al. 2012), and Binks et al. (2015) are quite similar, and compute proximity to the UVW distributions of stars using the formula

$$\begin{eqnarray}{\chi }^{2} & = & \displaystyle \frac{1}{3}\left(\displaystyle \frac{{({U}_{\mathrm{star}}-{U}_{\mathrm{NYMG}})}^{2}}{{\sigma }_{U,\mathrm{star}}^{2}+{\sigma }_{U,\mathrm{NYMG}}^{2}}\right.\\ & & +\displaystyle \frac{{({V}_{\mathrm{star}}-{V}_{\mathrm{NYMG}})}^{2}}{{\sigma }_{V,\mathrm{star}}^{2}+{\sigma }_{V,\mathrm{NYMG}}^{2}}\\ & & +\left.\displaystyle \frac{{({W}_{\mathrm{star}}-{W}_{\mathrm{NYMG}})}^{2}}{{\sigma }_{W,\mathrm{star}}^{2}+{\sigma }_{W,\mathrm{NYMG}}^{2}}\right),\end{eqnarray} \tag{ 11 }$$

which computes the difference between the star's UVW and the moving group's UVW in units roughly analogous to standard deviations. This technique does not deal with missing information as easily as a convergence-style routine; it is impossible to calculate the UVW without full kinematic information, requiring the use of estimates or repeated calculations over a range. Differences are mostly in implementation. Shkolnik et al. (2012) studied 14 NYMGs and assumed a 2 km s⁻¹ velocity dispersion for all groups; they required a maximum χ² value of six and a maximum velocity difference (to avoid identifying stars with large uncertainties) of 5 km s⁻¹. They calculated photometric distances where parallaxes were not available; all targets had γ. Binks et al. (2015) studied 10 groups with individualized dispersions and set a χ² cutoff of 3.78, and had a maximum velocity difference of 5 km s⁻¹. Riedel et al. (2014) studied 13 groups with individualized dispersions, and accepted a maximum χ² (there, called γ) value of 4.0. There was no limit on actual velocity difference. Where radial velocity did not exist, they computed UVW values for a range of radial velocities from −100 to +100 km s⁻¹ and took the minimum χ² (and corresponding RV) from a polynomial fit as the correct value; all targets had parallaxes.

The kinematic moving group identification technique in the SACY papers (e.g., Torres et al. 2008, p. 757) is similar to, but simpler than the Shkolnik et al. (2012), Riedel et al. (2014), and Binks et al. (2015) method, with $F$ $=\,(p\,\times {({M}_{{\rm{v}}}-{M}_{{\rm{v}},\mathrm{iso}})}^{2}$ $+\ {({U}_{\mathrm{star}}-{U}_{\mathrm{NYMG}})}^{2}$ $+\ {({V}_{\mathrm{star}}-{V}_{\mathrm{NYMG}})}^{2}$ $+{({W}_{\mathrm{star}}-{W}_{\mathrm{NYMG}})}^{2}{)}^{\tfrac{1}{2}}$, where ${M}_{{\rm{v}}}-{M}_{{\rm{v}},\mathrm{iso}}$ is the difference between the absolute V magnitude of the star and the absolute V magnitude of a star of a particular V − I_C color, on the age-appropriate isochrone, and F is the velocity-space separation between a star (utilizing a best-fit UVW for the star if the distance is unknown) and a moving group. In practice, they claim that the scaling constant p was only non-zero for the Octans moving group, and therefore the equation is generally a simple velocity separation independent of uncertainties. The cutoff of good values starts at F = 3.5, but the process used by SACY is to iteratively minimize UVW for a cluster of young stars until all outliers are removed. The team has considered nine groups: β Pic, Tuc-Hor, Columba, Carina, TW Hya, Cha, Octans, Argus, and AB Dor.

The code used by Klutsch et al. (2014) considers membership in two ways: first, a method using Gaussian representations of the UVW velocities (and their errors, formed by error propagation); second, member-to-member analysis, wherein the stars are compared to both the moving group center and the other known members; this allows for more variation and is potentially more robust at dealing with non-Gaussian distributions of stars.

We present the results of comparing LACEwING to the BANYAN, BANYAN II, and Convergence codes in Table 8. For each group, we record the number of members in our bona fide sample (Column 2), the number of objects identified by the code as members of that group (Column 3), and the number of objects identified as members that were genuine members (Column 4), with the corresponding number of false positives (Column 5).

Table 7.  Moving Group Kinematic Properties

		Values			Axes			Angles
Group	Number	U	V	W	A	B	C	UV	UW	VW
Name	BAFGK	(km s−1)	(km s−1)	(km s−1)	(km s−1)	(km s−1)	(km s−1)	(rad)	(rad)	(rad)
Cha	17	−10.9	−20.4	−9.9	0.8	1.3	1.4	0	0	0
η Cha	6	−10.2	−20.7	−11.2	0.2	0.1	0.1	0	0	0
TW Hya	7	−10.954	−18.036	−4.846	3.043	2.332	1.703	0.227	0.022	0.098
β Pic	34	−10.522	−15.964	−9.2	3.167	2.039	1.609	0.020	0.045	0.238
32 Ori	12	−11.8	−18.5	−8.9	0.4	0.4	0.3	0	0	0
Octans	22	−13.673	−4.8	−10.927	1.749	1.678	1.029	0.32	−0.52	0.241
Tuc-Hor	63	−9.802	−20.883	−1.023	4.01	2.883	1.458	−0.042	−0.588	0.568
Columba	52	−12.311	−21.681	−5.694	2.321	1.43	1.322	0.470	−0.142	0.329
Carina	22	−10.691	−22.582	−5.746	1.763	0.532	0.178	0.341	0.092	0.044
Argus	38	−22.133	−12.122	−4.324	1.992	1.755	0.774	−0.088	−0.026	0.002
AB Dor	86	−7.031	−27.241	−13.983	2.136	1.929	1.859	0.041	0.050	0.182
Car-Near	10	−27.020	−18.255	−3.021	3.044	1.819	1.147	0.023	0.149	−0.286
Coma Ber	104	−2.512	−5.417	−1.204	1.868	1.364	1.876	0.057	0.106	−0.202
Ursa Major	55	14.278	2.392	−8.987	2.64	0.594	0.407	−0.799	−0.766	0.500
χ01 For	14	−12.29	−20.95	−4.9	0.98	0.92	1.07	0	0	0
Hyades	260	−41.1	−19.2	−1.4	0.23	0.23	0.23	0	0	0
Field (Sirius)	5800	8	2	−7.25	12	6	9	−0.436	0	0
Field (Coma1)	4400	−10	−8	−7.25	9	6	7	−0.436	0	0
Field (Coma2)	2700	15	−18	−7.25	14	7	7	−0.436	0	0
Field (Hyades)	5800	−32	−17	−7.25	12	6	9	0.5	0	0
Field (Pleiades)	5800	−12	−24	−7.25	10	6	9	−0.436	0	0
Field (ζ Her)	800	−35	−48	−7.25	14	6	8	−0.436	0	0
Field	14800	−11.1	−25	−7.25	50	25	25	0	0	0
Group		X	Y	Z	D	E	F	XY	XZ	YZ
Name		(pc)	(pc)	(pc)	(pc)	(pc)	(pc)	(rad)	(rad)	(rad)
Cha		54	−92	−26	3	6	7	0	0	0
η Cha		33.4	−81	−34.9	0.4	1	0.4	0	0	0
TW Hya		16.816	−51.33	21.194	17.9	7.681	5.197	−0.847	0.034	0.183
β Pic		7.075	−3.509	−16.277	30.736	16.323	7.186	−0.06	0.043	−0.337
32 Ori		−89.634	−29.47	−24.34	3.4	3.4	3.4	0	0	0
Octans		15.913	−95.179	−63.138	64.92	20.831	13.888	0.059	−0.107	0.08
Tuc-Hor		5.477	−19.146	−35.177	22.83	13.179	6.713	−0.175	−0.043	0.287
Columba		−27.056	−26.369	−31.674	22.794	24.357	15.479	0.497	0.074	0.317
Carina		18.582	−65.598	−21.795	16.127	12.938	3.63	−0.747	0.27	−0.24
Argus		16.075	−21.027	−3.39	28.119	25.709	6.084	−0.434	−0.104	−0.018
AB Dor		−4.323	1.391	−17.372	23.857	21.531	15.014	−0.16	0.085	0.221
Car-Near		−2.961	−19.919	−1.955	23.641	10.619	5.037	−1.176	−0.532	1.006
Coma Ber		−6.706	−6.308	87.522	10	10	10	0	0	0
Ursa Major		−6.704	10.134	21.622	1.388	0.763	0.251	0.559	0.312	0.682
χ01 For		−28.3	−46.3	−83.4	2.9	2.9	2.9	0	0	0
Hyades		−43.1	0.7	−17.3	10	10	10	0	0	0
Field (Sirius)		−0.18	2.1	3.27	120	120	120	0	0	0
Field (Coma1)		−0.18	2.1	3.27	120	120	120	0	0	0
Field (Coma2)		−0.18	2.1	3.27	120	120	120	0	0	0
Field (Hyades)		−0.18	2.1	3.27	120	120	120	0	0	0
Field (Pleiades)		−0.18	2.1	3.27	120	120	120	0	0	0
Field (ζ Her)		−0.18	2.1	3.27	120	120	120	0	0	0
Field		−0.18	2.1	3.27	120	120	120	0	0	0

Download table as:  ASCIITypeset image

For instance, LACEwING (in young star mode) identifies 24 Cha members in the bona fide sample, and all 24 of its recovered members are known members of Cha: a 0% false positive rate, but not a 100% recovery rate given that there were 27 bona fide members of Cha to find. In contrast, LACEwING (in both modes) did recover all 11 bona fide members of χ⁰¹ For, and did not identify other stars as χ⁰¹ For members.

Looking at false positive rates is not as useful for the lithium sample (Table 9), because 518 out of the 930 lithium-rich stars are not previously known members of moving groups and many are likely to actually be members that have not been investigated kinematically before.

Table 8.  Comparison of Completion in the Bona Fide Data Sets

Group	Real	Identified	Recovered	False Positive	Code
(1)	(2)	(3)	(4)	(5)	(6)
Cha	27	24	24	0	L-Y
		22	22	0	L-F
η Cha	2	2	2	0	L-Y
		2	2	0	L-F
TW Hya	17	15	15	0	L-Y
		10	10	0	L-F
		16	16	0	B1
		16	16	0	B2-Y
		16	16	0	B2-F
		19	6	13	C
β Pic	28	22	22	0	L-Y
		14	14	0	L-F
		43	27	16	B1
		25	25	0	B2-Y
		25	25	0	B2-F
		38	14	24	C
32 Ori	10	9	9	0	L-Y
		8	8	0	L-F
Octans	45	44	44	0	L-Y
		17	17	0	L-F
Tuc-Hor	32	30	29	1	L-Y
		31	30	1	L-F
		30	29	1	B1
		31	30	1	B2-Y
		31	30	1	B2-F
		70	28	42	C
Columba	16	14	13	1	L-Y
		10	10	0	L-F
		29	15	14	B1
		15	14	1	B2-Y
		15	14	1	B2-F
		34	8	26	C
Carina	4	4	4	0	L-Y
		3	3	0	L-F
		4	4	0	B1
		7	4	3	B2-Y
		7	4	3	B2-F
Argus	6	8	5	3	L-Y
		5	5	0	L-F
		8	6	2	B1
		7	6	1	B2-Y
		7	6	1	B2-F
AB Dor	35	33	30	3	L-Y
		26	26	0	L-F
		35	35	0	B1
		30	30	0	B2-Y
		30	30	0	B2-F
		74	28	46	C
Car-Near	9	3	3	0	L-Y
		2	2	0	L-F
Coma Ber	45	33	33	0	L-Y
		31	31	0	L-F
Ursa Major	5	5	5	0	L-Y
		5	5	0	L-F
χ01 For	11	11	11	0	L-Y
		10	10	0	L-F
Hyades	0	0	0	0	L-Y
		0	0	0	L-F

Note. Comparison of the recovery rates of LACEwING in Young star (L-Y) and Field star (L-F) calibration to BANYAN (B1), BANYAN II in Young star (B2-Y) and Field star (B2-F) calibration, and the Convergence (C) code. In this context, "membership" means the highest probability was for the group in question.

Download table as:  ASCIITypeset images: 1 2

The comparisons in Tables 8 and 9 are not "recoveries" in the traditional sense, as there are no widely accepted correct answers (apart from, perhaps, the bona fide member list). LACEwING, BANYAN, and BANYAN II are roughly tied in terms of accuracy, except with the AB Dor moving group, where BANYAN clearly outperforms all other codes.

Tests of LACEwING while under development show that as more groups are added to the kinematic model, the recovery rate drops, but the false positive rate drops as well. While developing LACEwING, attempts were made to strike a compromise between recovery rates and false positive rates.

Table 9.  Comparison of Completion in the Lithium Data Sets

Group	Real	Identified	Recovered	False Positive	Code
(1)	(2)	(3)	(4)	(5)	(6)
Cha	30	30	25	5	L-Y
		27	24	3	L-F
η Cha	15	17	14	3	L-Y
		15	14	1	L-F
TW Hya	31	28	25	3	L-Y
		11	11	0	L-F
		34	27	7	B1
		28	27	1	B2-Y
		28	27	1	B2-F
		65	12	53	C
β Pic	44	24	19	5	L-Y
		13	12	1	L-F
		77	36	41	B1
		38	26	12	B2-Y
		37	26	11	B2-F
		78	20	58	C
32 Ori	0	1	0	1	L-Y
		0	0	0	L-F
Octans	22	31	21	10	L-Y
		9	9	0	L-F
Tuc-Hor	68	61	55	6	L-Y
		46	43	3	L-F
		54	52	2	B1
		53	50	3	B2-Y
		53	50	3	B2-F
		94	49	45	C
Columba	41	35	15	20	L-Y
		12	8	4	L-F
		70	32	38	B1
		33	19	14	B2-Y
		33	19	14	B2-F
		79	14	65	C
Carina	21	19	8	11	L-Y
		8	4	4	L-F
		17	8	9	B1
		19	7	12	B2-Y
		19	7	12	B2-F
Argus	28	23	7	16	L-Y
		8	4	4	L-F
		24	11	13	B1
		12	6	6	B2-Y
		11	6	5	B2-F
AB Dor	65	51	36	15	L-Y
		25	23	2	L-F
		62	50	12	B1
		27	24	3	B2-Y
		25	24	1	B2-F
		139	48	91	C
Car-Near	4	3	1	2	L-Y
		0	0	0	L-F
Coma Ber	0	1	0	1	L-Y
		0	0	0	L-F
Ursa Major	11	3	1	2	L-Y
		1	1	0	L-F
χ01 For	1	3	1	2	L-Y
		1	1	0	L-F
Hyades	0	11	0	11	L-Y
		0	0	0	L-F

Note. Same as Table 8, but for the lithium sample.

Download table as:  ASCIITypeset images: 1 2

This group was discovered by Mamajek et al. (2000) during examination of the η Cha open cluster. For a time, it was assumed (e.g., Torres et al. 2008, p. 757) to be co-eval with the open cluster and likely part of the same star-forming event; Murphy et al. (2013) recently used evolutionary models that suggest it is younger than η Cha. The hottest (and likely most massive) member of Cha is  Cha itself, a B9V star.

First discovered by Mamajek et al. (1999), η Cha is the smallest and third-closest open cluster to the Sun. The hottest member of the open cluster is η Cha itself, a B9V star whose companions form most of the members of the open cluster. Only two members of η Cha (η Cha and RS Cha) have parallaxes, so our moving group parameters are taken from Murphy et al. (2010).

Discovered by de la Reza et al. (1989), TW Hya was the second known NYMG discovered (after Ursa Major), and TW Hya itself is the closest classical T Tauri star to the Sun. Its UVW velocity and projected sky position are similar to the Lower Centaurus Crux region of the Sco-Cen star-forming complex, and the group itself extends from roughly 30 pc from the Sun to the near edge of Lower Centaurus Crux. Many proposed TW Hya members are now thought to actually be part of that background group. The hottest member of TW Hya is not actually TW Hya itself (K6Ve), but TWA 11 (HR 4796), an A0V star.

This group was first noticed by Mamajek (2007) as a small knot of stars around the B5V+B7V binary star 32 Ori; very little has been studied about this group of stars. Where most of the NYMGs have been linked to origins in the Sco-Cen star-forming region, 32 Ori may have a different origin: Bouy & Alves (2015) proposed a different arrangement of gas near the Sun where instead of Gould's belt, there are three parallel streams, one of which would stretch from 32 Ori to the Orion Nebula Complex.

Discovered by Barrado y Navascués et al. (1999b), β Pic is one of the moving groups that effectively surrounds the Sun. The most massive member of the moving group is HR 6070 (A0V), though our tracebacks actually rejected it as a member. The hottest member of β Pic not rejected by our tracebacks is β Pic itself, an A5V star with a planet and prominent debris disk.

Discovered by van den Ancker et al. (2000), Capricornus had only two proposed members, one of which (BD-17 6127 AB) is now identified as a β Pic member; the other, HD 356823, matches no known group. No attempt was made to consider this moving group for inclusion in LACEwING.

First published by Zuckerman & Song (2004), most of this moving group's members are now thought to be part of Argus and Cha. This group was never considered for inclusion in LACEwING.

First proposed by Torres et al. (2008, p. 757), Octans is one of the most distant moving groups. Based on our ellipse fits, it is also the largest (and therefore least dense) moving group, and LACEwING has difficulty differentiating its members from field stars. Very few papers have studied Octans, apart from Murphy & Lawson (2015). The hottest star in Octans is HD 36968, an F2 star. Every other group except Carina and Car-Near has a hotter member in the B5-A5 range, suggesting that there is either something different about Octans, or hotter members remain to be identified.

The Horologium moving group was discovered by Torres et al. (2000), followed by the Tucana moving group (Zuckerman et al. 2001), and then the realization that they were two parts of the same group (Zuckerman 2001). Thanks primarily to the work of Kraus et al. (2014a), Tuc-Hor has the most known members of any NYMG, though most are M dwarfs. Considering only the BAFGK members, Tuc-Hor is still likely smaller (63 BAFGK members) than AB Dor (86 BAFGK members). The hottest star in Tuc-Hor is α Pav, a B2IV star, although the β Tuc sextuple system (B9V+A0V+A2V+A7V+unknowns) may be the most massive.

For most NYMGs, their members have generally non-existent probabilities of membership in other NYMGs. This is not true for Tuc-Hor: a large fraction of Tuc-Hor members also have low but significant probabilities of membership in Columba, and somewhat less significant probability of membership in β Pic.

The Columba moving group was first announced by Torres et al. (2008, p. 757) along with the Carina moving group, as subdivisions of a larger "GAYA" complex that also contained Tuc-Hor. While Columba does have a different UVW velocity from Tuc-Hor and Carina, its spatial location at one end of the Tuc-Hor (see Figure 13) and persistently similar age to the other two groups continue to suggest that this may not be distinct from Tuc-Hor. LACEwING finds that many Columba bona fide members have low probabilities of membership in Tuc-Hor. The most massive member of Columba is HR 1621 (B9V). The known planet host (Marois et al. 2008) and putative Columba member (Baines et al. 2012) HR 8799 survived both TRACEwING filtering and LACEwING analysis as a bona fide member of Columba, despite the fact that its high mass (as an A5V star) and position far from other known members are at odds with the idea of mass segregation.

As mentioned in Section 7.10, Carina was first discovered by Torres et al. (2008, p. 757) as part of the "GAYA" complex, a larger group that included Tuc-Hor and Columba; like Columba, Carina sits physically adjacent to Tuc-Hor. A preliminary version of LACEwING used in Riedel (2016b) and Faherty et al. (2016) actually excluded Carina because, with an earlier version of the Catalog, tracebacks removed all but two stars from Carina (see Section 7.17). Carina has a spatial volume only slightly larger than η Cha, Cha, and Ursa Major. The most massive member of Carina is HD 83096, an F0V star (which has a lower mass than most moving groups' largest members; see also Section 7.8).

Most members of Carina have low, but non-zero, probabilities of membership in Columba or Tuc-Hor. It is also a fairly poorly recovered group by all moving group codes studied; the best recovery in the lithium sample was 8 of 21 members, by both BANYAN and LACEwING in Young Star mode.

Discovered by Makarov & Urban (2000), Car-Vel was suspected of being related to the IC 2391 open cluster in much the same way as the IC 2391 Supercluster (Eggen 1991), although they share no members. Most of the stars thought to be in Car-Vel are now assigned as field stars, members of Carina, or members of Argus. This group was never considered for inclusion in LACEwING.

Argus was first identified by Torres et al. (2008, p. 757) as an improvement on the Car-Vel moving group, with which it shares some members. Argus shares a UVW velocity with the IC 2391 open cluster and is thought to be the product of the same star-forming event (De Silva et al. 2013). This association is problematic: with the cluster at a distance of ∼150 pc (Caballero & Dinis 2008), it is not clear how the stars could have reached the vicinity of the Sun yet have space velocities parallel to IC 2391 at the present day. Argus is also problematic in that Bell et al. (2015) found the members to not be co-eval and failed to compute an age for the NYMG. The hottest member of Argus is  Pav (A0V).

The problematic nature of Argus may also be exhibited in the lithium sample study, where none of the codes managed to find even half of the 28 lithium-detected members of Argus.

AB Dor was first identified by Torres et al. (2003) as "AnA" and by Zuckerman et al. (2004) as AB Dor (and possibly pre-discovered by Asiain et al. 1999 as subgroup B4.) AB Dor has a space velocity distribution encompassing that of the Pleiades, and has been thought to be a product of the same star-forming event (Luhman et al. 2005; Ortega et al. 2007). Though initially thought to be a younger population at 50 Myr (Close et al. 2005), it is now believed to be as old as 150 Myr (Bell et al. 2015), which would make it older than the Pleiades. AB Dor is the largest of the NYMGs based on the number of systems with BAFGK primaries. However, AB Dor is old enough that lithium and surface gravity are less effective indicators for low-mass objects, and chemical tagging studies (Barenfeld et al. 2013; McCarthy & Wilhelm 2014) have determined that perhaps as many as half of the AB Dor members may be interlopers. AB Dor is effectively an all-sky moving group. The hottest member of AB Dor is Alnair (B6V).

LACEwING recovers all but eight bona fide members of AB Dor in Field star mode and all but five members in young star mode; in all cases the rejected stars match no other moving group. In both cases, it failed to recover Alnair and δ Scl (the hottest member that it does recover is HR 1014, A3V).

Car-Near was identified by Zuckerman et al. (2006) as a nearby older population of stars. The hottest member of Car-Near is HR 3070 (F1), which is cooler than most group's hottest members (see Sections 7.8 and 7.11).

Car-Near is not well recovered by LACEwING due to its large volume and small membership. One member (LP 356-14) is identified in both young and field star modes as a member of Argus, the group whose UVW velocity Car-Near most closely resembles. HR 3070 is an Argus member in Young Star mode. The rest match either Car-Near or no group at all. Despite this poor recovery, Car-Near still appears to be a real group. It does not have an unusually large present-day volume like Oct-Near (Section 7.16) and a corresponding complete failure of recovery, its members are largely only members of the group itself, and it does produce self-consistent tracebacks, unlike Her-Lyr (Section 7.17). It seems most likely that other members of Car-Near will be found, increasing the spatial density and therefore the recovery rate.

Oct-Near was identified by Zuckerman et al. (2013) as a potential very-nearby association of stars (including EQ Peg AB at 6.2 pc) with similar velocities to Octans and Castor. Though considered for inclusion in LACEwING, Oct-Near posed two problems for inclusion. First, as presented in Zuckerman et al. (2013), the group had multiple apparent ages between 30 and 200 Myr. Second, the present-day moving group ellipses were more than three times the size of the other moving groups and the resulting groups so sparse that LACEwING could not recover any of the supposed bona fide members. Accordingly, it is not included in this implementation of LACEwING. The hottest member of Octans-Near was 34 Psc (B9 V).

Her-Lyr was first identified by Gaidos (1998) and Fuhrmann (2004), and was comprised almost entirely of nearby stars; it would have been another all-sky moving group. The existence of Her-Lyr has been disputed for some time; as Mamajek (2015) notes, multiple papers (most recently, Eisenbeiss et al. 2013) have identified members of Her-Lyr, but none have consistent lists of members. Nevertheless, in recognition of the fact that the supposed members of Her-Lyr were a population of lithium-rich stars older than AB Dor, we attempted to retain Her-Lyr. Her-Lyr did appear in the preliminary LACEwING calibration used in Riedel (2016b) and Faherty et al. (2016), but with data from our updated Catalog of Suspected Nearby Young Stars (Section 4), only one star—V0439 And—remained within the boundaries of Her-Lyr and survived traceback filtering. Her-Lyr is not included in this implementation of LACEwING. The hottest member of Her-Lyr was α Cir (A7V), which was not listed as a bona fide member in Eisenbeiss et al. (2013).

Castor was first identified by Anosova & Orlov (1991) in a relatively wide (6 km s⁻¹) search for potential clusters around prominent triple star systems. Though refined in multiple papers over the ensuing 25 years, the existence of Castor has been questioned and debunked by Mamajek et al. (2013) on the grounds that the most prominent members (Vega, Fomalhaut, LP 944-020, and Castor) were nowhere near each other even 10 Myr ago and could not have formed in the same molecular cloud. This determination was made with a different data set and a simple linear traceback. The existence of Castor was also disputed by Zuckerman et al. (2013) on the basis of differing ages and large velocity spreads within the group.

We have attempted to trace back Castor ourselves, starting with all 84 members ever identified as Castor members (regardless of current identification) from the Catalog of Suspected Nearby Young Stars with parallaxes and radial velocities. Only 33 stars survived three rounds of TRACEwING filtering, at which point the group included neither Castor itself nor LP 944-020, and was still the largest moving group in terms of volume at formation (Figure 17). Castor's volume at formation is indistinguishable from the size of the fake moving group of field stars. Castor, a sextuple system with four A-type stars, is likely both the hottest and most massive member of the moving group.

We cannot conclusively rule out the existence of Castor because it is still smaller than our simulated moving group between the ages of 200–400 Myr (covering both the age given by Barrado y Navascues 1998 and the ages of the prominent members collected in Mamajek et al. 2013), but we suspect this is further evidence that Castor is not a real moving group, and we have not included it in LACEwING.

The Ursa Major moving group, often referred to as the Sirius Supercluster in older literature and occasionally referred to as a cluster (Mamajek 2015), was the first moving group discovered, by R. A. Proctor in the late 19th century. Castro et al. (1999) published a chemical tagging analysis that revealed that the group members are identifiably rich in barium. King et al. (2003) published the most recent large-scale study of the group and determined that Sirius is not a likely member. Ursa Major is represented here by only the core members from King et al. (2003) and is thus one of the smallest moving groups in LACEwING. The hottest member of Ursa Major is  UMa (A0pCr), although Mizar-Alcor, a sextuplet with five A-type stars, is likely to be more massive.

Ursa Major is another troubling group; the final traceback set of Ursa Major members did not include Mizar-Alcor, and the resulting LACEwING calibration missed 8 of the 11 lithium-rich members in young mode (although LACEwING does recover Mizar-Alcor as a member).

Coma Ber (Melotte 111, Collinder 256) is an open cluster 86 pc distant consisting of roughly 195 stars from our limited literature search, of which 104 are BAFGK members. Thus, Coma Ber is second only to the Hyades in terms of size. Coma Ber was a difficult moving group to add to our simulation of the Solar Neighborhood because it has a very low UVW velocity and the simulation produced proper motions indistinguishable from zero for many of the simulated members. Fortunately, LACEwING's field population prevents LACEwING from identifying large portions of the sky as Coma Ber members. The hottest member of Coma Ber is AI Com (A0p).

χ⁰¹ For (also known as Alessi 13) was first published in a catalog by Dias et al. (2002). It has remained obscure for the past decade, but appears to be real (E. E. Mamajek 2016, private communication). The only available age estimates are from Pöhnl & Paunzen (2010) and Kharchenko et al. (2013), which hover around 500 Myr, though Mamajek believes the group may be younger (Mamajek 2016). For the purposes of traceback we intended to use 450–550 Myr, but only three members of the group have parallaxes. We used the UVW properties from the private communication with E. E. Mamajek instead, and calculated the XYZ positions from the α, δ, Dist., and tidal radius. The hottest member of χ⁰¹ For is χ⁰¹ For (A1V) itself.

The Hyades open cluster (Melotte 25, Collinder 50) has been known since antiquity. Based on analyses of Röser et al. (2011) and Goldman et al. (2013), there are 724 known members of the Hyades, of which 260 are BAFGK members (including giants and white dwarfs). The hottest stellar member (excluding white dwarfs) is θ⁰² Tau, an A7III giant. The age is believed to be between 600 (Zuckerman & Song 2004) and 800 Myr (Brandt & Huang 2015), and is the upper limit of what we consider to be a young group.

The Hyades is represented here by a conversion of the properties in van Leeuwen (2009) to UVWXYZ.

A large fraction of our lithium sample—582 of the 930 star systems—do not trace back to any of the NYMGs in young star mode. Twenty-four of these just miss (>10% probability of membership in at least one group) the usual 20% probability cut, which could be due to only having proper motion available for membership assignment, but that still leaves 534 of 930 systems (57%) as nonmembers of the groups. This behavior is not unique to LACEwING. Running the same sample through BANYAN, which had the highest recovery rate of any of the codes, 589 stars fell into the "Old" category. BANYAN II in field star mode identified 369 stars as young field objects and 355 stars as old field objects; in young star mode it identified 720 members of the lithium sample as "Old." This is not a perfect comparison, as the BANYAN codes test for fewer groups than LACEwING, and some of the lithium-detected stars are members of groups older than BANYAN's oldest.

If we break down the entire lithium sample, only 487 of the star systems have ever been considered as potential members of groups that LACEwING tests for. An additional 79 have only been considered as members of moving groups LACEwING does not test for (Her-Lyr, Oct-Near, Castor, Car-Vel, Local Association, Hyades Supercluster, IC 2391 Supercluster). The remaining 364 systems have never been considered as members of any specific young group.

This general behavior has been noted before, particularly in surveys that did not filter by kinematic match before following up stars for further observations (Shkolnik et al. 2009, 2012; Riedel et al. 2014, 2016b; Binks et al. 2015). These results have suggested that there is an unidentifiable population of young stars.

If we accept that these stars are young, there are five possibilities for their origins.

• 1.  
  The stars did form as part of the known NYMGs, and flaws in our data or models are responsible for the stars that are not associated with a group.
• 2.  
  The stars did form as part of the known NYMGs, but dynamical interactions (most likely an ejection from a higher-order multiple star system early in its evolution) gave them high relative velocities such that they do not kinematically match the group that they formed with. Chemical tagging would be useful for identifying such systems, provided it is possible to uniquely identify a moving group in that way. Kinematic tracebacks that place a star near another member of a moving group of the appropriate age would be suggestive, but would require advancements in both data precision and technique from what is presented here.
• 3.  
  The stars did form as part of groups that are known but not nearby. LACEwING only tests for groups currently known to extend within 100 pc of the Sun. The Scorpius–Centaurus and Taurus–Auriga star-forming regions are less than 200 pc away. With a typical velocity dispersion of 1.5 km s⁻¹ (Preibisch & Mamajek 2008), a star could move from 118 pc (the canonical distance to the ∼16 Myr old Lower Centaurus Crux region of Scorpius–Centaurus; Preibisch & Mamajek 2008) to within 100 pc of the Sun in just over 10 Myr.
• 4.  
  The stars did form as part of unknown NYMGs that we have not yet discovered. 32 Ori and χ⁰¹ For are relatively unexplored groups; the All Sky Young Association (Torres et al. 2016) has been announced but no particulars have yet been given. Other, smaller groups that have no members more massive than M dwarfs may yet be hiding in the Solar Neighborhood.
• 5.  
  The stars did not form as part of any groups at all, and were instead part of a one-off star formation event. The NYMGs already range greatly in size from the η Cha open cluster (21 known systems as of 2015 January) to the Tuc-Hor moving group (209 known systems as of 2015 January); it is not clear what the smallest star formation event can be.

The LACEwING algorithm is available as a function that can be called from other programs (see Appendix A.2), but if it is run directly from the command line it defaults to reading from a comma-separated value file:

$$\begin{eqnarray*}&&{\mathtt{python}}\ {\mathtt{lacewing}}.{\mathtt{py}}\ {\mathtt{inputfile}}\ [{\mathtt{calibration}}]\\ &&[{\mathtt{output}}\ {\mathtt{filename}}]\ [{\mathtt{verbose}}\ {\mathtt{output}}]\ [{\mathtt{g}}.{\mathtt{o}}.{\mathtt{f}}]\end{eqnarray*}$$

LACEwING (using the astropy.io.ascii module) looks for a one-line header with any or all of the columns selected below (broadly, they are either common names for the quantities or AAS Journal standards). These can be given in any order; any other headers present in the file are ignored. Most of the errors are optional and if they (and their header) are not present, LACEwING will substitute default values.

• 1.  
  Name—ascii, any length; should not contain commas itself.
• 2.  
  RA, DEC (or RAdeg, DEdeg)—right ascension and declination in decimal degrees, J2000/ICRS equinox.
• 3.  
  RAh, RAm, RAs, DE-, DEd, DEm, DEs—sexagesimal coordinates, J2000/ICRS equinox, split into seven columns with a separate declination +/− flag. These will only be read if RA and DEC are invalid, empty, or do not exist.
• 4.  
  eRA, eDEC (or e_RAdeg, e_DEdeg)—${\rm{R}}.{\rm{A}}.\mathrm{cosdecl}.$ and decl. uncertainties in milliarcseconds. These, and their header keyword, are optional and will default to 1000 milliarcseconds.
• 5.  
  pmRA, pmDEC (or pmRA, pmDE; or pmra and pmdec)—${\mu }_{{\rm{R}}.{\rm{A}}.\cos \mathrm{decl}.}$ and μ_decl. in mas yr⁻¹, J2000/ICRS equinox.
• 6.  
  epmRA, epmDEC (or e_pmRA,e_pmDE; or epmra, epmdec)—$e{\mu }_{{\rm{R}}.{\rm{A}}.\cos \mathrm{decl}.}$ and eμ_decl. uncertainties in mas yr⁻¹, J2000/ICRS equinox. These, and their header, are optional and will default to 10 mas yr⁻¹.
• 7.  
  pi (or plx)—trigonometric parallax in milliarcseconds. It is highly recommended that the photometric parallax estimates are not used.
• 8.  
  epi (or e_plx)—trigonometric parallax uncertainty in milliarcseconds. This must be present in order to use the parallax.
• 9.  
  rv (or HRV)—radial velocity in km s⁻¹.
• 10.  
  erv (or e_HRV)—radial velocity uncertainty in km s⁻¹. This must be present in order to use the radial velocity.
• 11.  
  Note—any text you wish to have duplicated in the file output, such as a previously known membership assignment.

The additional command-line parameters are

• 1.  
  calibration—type "young" to switch to the alternative LACEwING calibration where stars are assumed to be young. Defaults to "field"
• 2.  
  output filename—enter an output filename, otherwise LACEwING defaults to lacewing_output.csv.
• 3.  
  verbose—If this is "verbose," LACEwING will output a more detailed report as described below. It defaults to a much more compact output format.
• 4.  
  g.o.f—If this is set to anything other than "percentage," LACEwING outputs the combined goodness-of-fit statistic rather than determining the moving group membership probability.

The LACEwING output format, in the regular version, is a comma-separated value file containing, in order:

• 1.  
  Star name
• 2.  
  Note (directly from the input)
• 3.  
  Best-matching moving group
• 4.  
  Membership probability for the best-matching moving group (unless the "g.o.f" flag is set on the command line)
• 5.  
  Kinematic distance (in parsecs) for the best-matching group
• 6.   
  Kinematic distance uncertainty (in parsecs) for the best-matching group
• 7.   
  Kinematic radial velocity (in km s⁻¹) for the best-matching group
• 8.  
  Kinematic radial velocity uncertainty (in km s⁻¹) for the best-matching group
• 9.   
  Membership probability for Cha
• 10.  
  Membership probability for η Cha
• 11.  
  Membership probability for TW Hya
• 12.  
  Membership probability for β Pic
• 13.  
  Membership probability for 32 Ori
• 14.  
  Membership probability for Octans
• 15.  
  Membership probability for Tuc-Hor
• 16.  
  Membership probability for Columba
• 17.  
  Membership probability for Carina
• 18.  
  Membership probability for Argus
• 19.   
  Membership probability for AB Dor
• 20.  
  Membership probability for Carina-Near
• 21.  
  Membership probability for Coma Ber
• 22.   
  Membership probability for Ursa Major
• 23.  
  Membership probability for χ⁰¹ For
• 24.  
  Membership probability for Hyades

The Verbose Output formula is drastically different, consisting of a block of 16 lines detailing the match between the star and every group, one at a time. It is most useful for tracking down the reason why a particular star did not match a particular group.

• 1.  
  Name
• 2.  
  R.A. (degrees, from the input)
• 3.  
  Decl. (degrees, from the input)
• 4.  
  Moving group
• 5.  
  Text string: "PROB=" (or "SIG=," if the g.o.f. flag has been set on the command line)
• 6.  
  Membership probability in the moving group (or combined goodness of fit if the g.o.f. flag has been set)
• 7.  
  Text string: "PM="
• 8.  
  Proper motion goodness-of-fit metric
• 9.  
  Kinematic (predicted) ${\mu }_{{\rm{R}}.{\rm{A}}.\cos \mathrm{decl}.}$ (mas yr⁻¹)
• 10.  
  Kinematic ${\mu }_{{\rm{R}}.{\rm{A}}.\cos \mathrm{decl}.}$ uncertainty (mas yr⁻¹)
• 11.  
  Kinematic μ_decl. (mas yr⁻¹)
• 12.  
  Kinematic μ_decl. uncertainty (mas yr⁻¹)
• 13.  
  Measured ${\mu }_{{\rm{R}}.{\rm{A}}.\cos \mathrm{decl}.}$ (mas yr⁻¹)
• 14.  
  Measured ${\mu }_{{\rm{R}}.{\rm{A}}.\cos \mathrm{decl}.}$ uncertainty (mas yr⁻¹)
• 15.  
  Measured μ_decl. (mas yr⁻¹)
• 16.  
  Measured μ_decl. uncertainty (mas yr⁻¹)
• 17.  
  Text string: "DIST="
• 18.  
  Distance goodness-of-fit metric
• 19.  
  Kinematic distance (parsecs)
• 20.  
  Kinematic distance uncertainty (parsecs)
• 21.   
  Measured distance ($\tfrac{1}{\pi }$, parsecs)
• 22.  
  Measured distance uncertainty ($\tfrac{{\sigma }_{\pi }}{{\pi }^{2}}$, parsecs)
• 23.   
  Text string: "RV="
• 24.  
  Radial velocity goodness-of-fit metric (km s⁻¹)
• 25.  
  Kinematic radial velocity (km s⁻¹)
• 26.  
  Kinematic radial velocity uncertainty (km s⁻¹)
• 27.   
  Measured radial velocity (km s⁻¹)
• 28.  
  Measured radial velocity uncertainty (km s⁻¹)
• 29.  
  Text string: "POS=" (or "KPOS=," if the spatial position is based on the kinematic distance)
• 30.  
  Spatial position goodness-of-fit metric
• 31.   
  3D separation between star and center of moving group (parsecs)
• 32.  
  3D separation uncertainty (parsecs)
• 33.  
  Note (directly from the input)

The lacewing_summary.py function will summarize the verbose output into the regular format.

$$\begin{eqnarray*}&&{\mathtt{python}}\ {\mathtt{lacewing}}\_{\mathtt{summary}}.{\mathtt{py}}\ {\mathtt{inputfile}}\end{eqnarray*}$$

The output file will be named inputfile.summary.csv

LACEwING as a function requires two calls. One, lacewing.moving_group_loader() returns a list of moving group classes loaded from the Moving_Group_all.csv file. This must be done before lacewing.lacewing(), so the resulting list can be passed to lacewing.lacewing(). This is done to avoid the redundancy of reloading the parameters every time.

lacewing.lacewing() accepts data on a single object at a time, and cannot be fed lists, arrays, or tuples. The returned list of moving group classes must be passed to lacewing.lacewing() every time. The first argument is required:

• 1.  
  moving_group—list, moving group class list (as created by lacewing.moving_group_loader()).

The remainder of the items are optional; if not present (or explicitly set to NoneType None) lacewing() will treat them as unknowns.

• 2.   
  young = string, if this is "young" then lacewing() will use the young star calibration. Otherwise (or if None, or not specified), lacewing() will use the field star calibration
• 3.   
  ra = float, decimal degrees, J2000/ICRS coordinates
• 4.  
  era = float, R.A. $\cos \,\mathrm{decl}.$ uncertainty in decimal degrees
• 5.  
  dec = float, decimal degrees, J2000/ICRS coordinates
• 6.  
  edec = float, decl. uncertainty in decimal degrees
• 7.   
  pmra = float, ${\mu }_{{\rm{R}}.{\rm{A}}.\cos \mathrm{decl}.}$ in arcsec yr⁻¹
• 8.  
  epmra = float, ${\mu }_{{\rm{R}}.{\rm{A}}.\cos \mathrm{decl}.}$ uncertainty in arcsec yr⁻¹
• 9.  
  pmdec = float, μ_decl. in arcsec yr⁻¹
• 10.   
  epmdec = float, μ_decl. uncertainty in arcsec yr⁻¹
• 11.   
  plx = float, π in arcseconds
• 12.  
  eplx = float, π uncertainty in arcseconds
• 13.  
  rv = float, RV in km s⁻¹
• 14.   
  erv = float, RV uncertainty in km s⁻¹

Note that lacewing() does not check to see if the uncertainties exist (or are not None), so (for example) supplying rv without erv will cause lacewing() to crash with a TypeError.

The output from lacewing() is a list of dicts, one dict per moving group in the order they appear in Moving_Group_all.csv. Each dict contains the following keys (which will be None unless data exists):

• 1.   
  "group"—string, group name
• 2.   
  "gof"—float, goodness-of-fit parameter
• 3.  
  "probability" float, membership probability (percent)
• 4.  
  "pmsig"—float, proper motion match metric
• 5.  
  "kin_pmra"—float, estimated ${\mu }_{{\rm{R}}.{\rm{A}}.\cos \mathrm{decl}.}$ (arcsec yr⁻¹)
• 6.  
  "kin_epmra"—float, estimated ${\mu }_{{\rm{R}}.{\rm{A}}.\cos \mathrm{decl}.}$ uncertainty (arcsec yr⁻¹)
• 7.   
  "kin_pmdec"—float, estimated μ_decl. (arcsec yr⁻¹)
• 8.   
  "kin_epmdec"—float, estimated μ_decl. uncertainty (arcsec yr⁻¹)
• 9.  
  "distsig"—float, distance match metric
• 10.  
  "kin_dist"—float, expected distance (parsecs)
• 11.   
  "kin_edist"—float, expected distance uncertainty (parsecs)
• 12.   
  "rvsig"—float, radial velocity match metric
• 13.   
  "kin_rv"—float, expected radial velocity (km s⁻¹)
• 14.   
  "kin_erv"—float, expected radial velocity uncertainty (km s⁻¹)
• 15.  
  "possig"—float, position match metric (using measured distance)
• 16.   
  "pos_esig"—float, position match metric uncertainty
• 17.   
  "pos_sep"—float, distance from moving group center (parsecs)
• 18.   
  "posksig"—float, position match metric (using kinematic distance)
• 19.  
  "pos_eksig"—float, kinematic position match metric uncertainty
• 20.   
  "pos_ksep"—float, kinematic distance from moving group center (parsecs)

There are two example implementations of lacewing.lacewing() in the repository: The default reader within the lacewing.py file, and the sample star generator in lacewing_montecarlo.py.

Step 1: Add the group to Moving_Groups_all.csv.

Generate new moving group parameters—Create a comma-separated value file in the format described above for the default csv loader, containing data on every member of the group.

Run this file through the moving group maker:

$$\begin{eqnarray*}&&{\mathtt{python}}\ {\mathtt{lacewing}}\_{\mathtt{mgpmaker}}.{\mathtt{py}}\ {\mathtt{inputfile}}\end{eqnarray*}$$

Also, run this through the 2D version to get parameters for lacewing_uvwxyz.py

$$\begin{eqnarray*}&&{\mathtt{python}}\ {\mathtt{lacewing}}\_{\mathtt{mgpmaker}}{\mathtt{2}}{\mathtt{d}}.{\mathtt{py}}\ {\mathtt{inputfile}}\end{eqnarray*}$$

The output files are .csv files named "Moving_Group_Group Name.dat" and "Moving_Group_Group Name_2d.dat" with U, V, W, A, B, C, UV, UW, VW, X, Y, Z, D, E, F, XY, XZ, and YZ values, with uncertainties (which LACEwING ignores) for each value on the second line. The first line is suitable for entering into Moving_Groups_all.csv. The 2D versions of those values should be stored in the A2, B2, C2, UV2 (etc.) columns. Note that UVW and XYZ values are the same for 2D and 3D projections, so there are no special 2D entries for them.

Alternatively, if you have the more customary UVW and XYZ values (no rotations) from a different source, you can enter those into the file, specifying 0 for all the rotation angles. The same values will apply for 3D and 2D.

Column "Name" should contain the name of the group, preferably less than 20 characters, with underscores instead of spaces.

Column "Number" is a best estimate of the current known members of the group.

Column "Weightednumber" is a cumulative fraction of stars in each group,

$$\begin{eqnarray}&&\displaystyle \frac{{\displaystyle \sum }_{0}^{i}\,{\mathrm{Number}}_{\mathrm{group}({\rm{i}})}}{{\displaystyle \sum }_{0}^{N}\,{\mathrm{Number}}_{\mathrm{group}({\rm{i}})}},\end{eqnarray} \tag{ 12 }$$

such that all are between 0 and 1, and should be generated from the number of members.

Column "uniform" specifies whether the group should be simulated as (0) a uniform spatial and Gaussian velocity distribution (unused), (1) a Gaussian spatial and Gaussian velocity distribution, like an open cluster, or (2) a spatial distribution with a scale height of 300 pc and a Gaussian velocity distribution, like the field population.

You can also add an Age (unused), References (unused), and the fractional RGB color components (for plots, such as lacewing_uvwxyz.py). Also, fill the membership coefficients section with zeros so that the program will function initially.

If you are removing a group, just delete its line and re-compute the "Weightednumber" column.

Step 2: Generate a new simulation—Generate a new Monte Carlo simulation.

$$\begin{eqnarray*}&&{\mathtt{python}}\ {\mathtt{lacewing}}\_{\mathtt{montecarlo}}.{\mathtt{py}}\ {\mathtt{iterations}}\\ &&{\mathtt{number}}\end{eqnarray*}$$

where "iterations" is the number of stars to draw out of the distributions now listed in Moving_Groups_all.csv. "Number" is a value to append to the names of the output files so that they will be kept separate, and is a way of more efficiently using computer resources.

Good values for "number" are the number of cores you are willing to devote to this process. Good values for "iterations" depends on the number of members in the least-populated group. This process needs roughly 1000 simulated members in total to populate the eventual membership histograms well enough to fit them. So, if you have a quad-core machine and the least-populated group accounts for fractionally 0.0001 of the stars, you need 10 million points and should generate 5 million with four running instances. An Intel Core i7 4700MQ running eight instances can generate about two million stars per day, and 8,000,000 entries takes up roughly 610 MiB of space per moving group.

The program will use random number generators to follow the procedure in Section 2.2 and generate simulated stars to run through lacewing.lacewing().

The output will be number files for each moving group, each containing iterations entries, with the goodness-of-fit values matching each star to the moving group using all seven possible combinations of input data, and a record of what group the star was generated as a member of.

Concatenate the files for each group into one file per group. It is preferable to move them to another subfolder.

Step 3: Generate and fit probability histograms—Run lacewing_percentages.py on the Monte Carlo output.

$$\begin{eqnarray*}&&{\mathtt{python}}\ {\mathtt{lacewing}}\_{\mathtt{percentages}}.{\mathtt{py}}\\ &&{\mathtt{montecarlofile}}\ [{\mathtt{young}}]\end{eqnarray*}$$

This program generates the histogram of "percentage of objects at each goodness-of-fit value that are actual members." Run it on one of the files from the Monte Carlo output to generate a wide variety of plots (as appearing in Section 2.2) and a file called group.percentages. This file contains the cumulative Gaussian fit parameters that need to be saved in Moving_Groups_all.csv. All of these files are stored in a folder called montecarlo, which is created if it does not exist.

If you specify "young," the program will generate the young star calibration; as it combs through the input file it will ignore all but an equal number of field stars to provide the case where the stars are assumed to be young. These calibrations (in group.young.percentages) should also be saved in Moving_Groups_all.csv, replacing the zeroes put there earlier.

You can now use LACEwING to predict members of your new moving group.

Requirements:

• 1.   
  A Python 2.7 interpreter with numpy, astropy, and matplotlib
• 2.   
  lacewing.py—the main routines necessary for LACEwING, specifically the csv loader
• 3.   
  lacewing_uvwxyz.py—the UVWXYZ fitting and plotting tool
• 4.  
  kinematics.py—routines to convert between equatorial (R.A., decl., π) and Galactic (XYZ) coordinates (and also kinematic tracebacks)
• 5.  
  ellipse.py—ellipse fitting and rotation routines
• 6.  
  astrometry.py—proper motion conversion routines and other important miscellaneous routines
• 7.  
  Moving_Groups_all.csv—a comma-separated value format file with the precalculated moving group parameters

A.4.1. General Usage

lacewing_uvwxyz.py calculates XYZ (where parallaxes exist) and UVWXYZ (where all six elements of kinematics exist) values, which are output to an output filename. It also generates 2D projection plots of stars plotted on the UVW and XYZ spatial axes, as shown in Figure 18:

$$\begin{eqnarray*}&&{\mathtt{python}}\ {\mathtt{lacewing}}\_{\mathtt{uvwxyz}}.{\mathtt{py}}\ {\mathtt{inputfile}}\ [{\mathtt{output}}\\ &&{\mathtt{filename}}]\ [{\text{}}{\mathtt{XYZ}}]\end{eqnarray*}$$

Input file format is the same as lacewing; it uses the same .csv loader).

The next two arguments are optional:

• 1.  
  output filename—enter an output filename, otherwise LACEwING defaults to lacewing_output.csv
• 2.  
  XYZ—if this string is "XYZ," the output .png image will have a second row of projected XYZ positions

There are two outputs:

• 1.  
  One image per star (named starname_decimal coordinates.png) with panels showing the U versus V, U versus W, and V versus W projections of 3D velocity (and, if "XYZ" was specified, X versus Y, X versus Z, and Y versus Z projections of 3D space) similar to Figure 18
• 2.  
  One output comma-separated value file with rows of Name, U, eU, V, eV, W, eW, X, eX, Y, eY, Z, and eZ values for each star.

• 1.  
  A Python 2.7 interpreter with numpy, astropy, and matplotlib
• 2.   
  tracewing.py—the kinematic traceback tool
• 3.   
  lacewing.py—for the csv loader
• 4.  
  kinematics.py—routines to trace the star back in time with epicyclic tracebacks
• 5.  
  ellipse.py—ellipse fitting and rotation routines
• 6.   
  astrometry.py—proper motion conversion routines and other important miscellaneous routines
• 7.  
  Moving_Group_Group Name_epicyclic.dat—stored parameter files for the moving groups. The files on GitHub have been calculated back to −800 Myr.

tracewing.py computes the traceback of stars to a given moving group, which must be present in a saved parameter file. The outputs are .png figures of the star traced back to the NYMG (Figure 7), covering the time between 0 and the end of the plotting time period.