January 18th 2023

Addressing Systematics in the Traceback Age of the
Systematics in the Traceback Age of the
Couture et al.

0000-0003-2604-3255]Dominic Couture
Planétarium Rio Tinto Alcan, Espace pour la Vie, 4801 av. Pierre-de Coubertin, Montréal, Québec, Canada
Trottier Institute for Research on Exoplanets, Université de Montréal, Département de Physique, C.P. 6128 Succ. Centre-ville, Montréal, QC H3C 3J7, Canada
Département de Physique, Université de Montréal, C.P. 6128 Succ. Centre-ville, Montréal, QC H3C 3J7, Canada

0000-0002-2592-9612]Jonathan Gagné
Planétarium Rio Tinto Alcan, Espace pour la Vie, 4801 av. Pierre-de Coubertin, Montréal, Québec, Canada
Trottier Institute for Research on Exoplanets, Université de Montréal, Département de Physique, C.P. 6128 Succ. Centre-ville, Montréal, QC H3C 3J7, Canada

0000-0001-5485-4675]René Doyon
Trottier Institute for Research on Exoplanets, Université de Montréal, Département de Physique, C.P. 6128 Succ. Centre-ville, Montréal, QC H3C 3J7, Canada
Département de Physique, Université de Montréal, C.P. 6128 Succ. Centre-ville, Montréal, QC H3C 3J7, Canada

We characterize the impact of several sources of systematic errors on the computation of the traceback age of the  (). We find that uncorrected gravitational redshift and convective blueshift bias absolute radial velocity measurements by , which leads to erroneously younger traceback ages by . Random errors on parallax, proper motion, and radial velocity measurements lead to an additional bias of  on traceback ages. Contamination of astrometric and kinematic data by kinematic outliers and unresolved multiple systems in the full input sample of  members and candidates of  also erroneously lowers traceback ages by . We apply our new numerical traceback analysis tool to a core sample of
carefully vetted members of  using  Data Release 3 (DR3) data products and other kinematic surveys. Our method yields a corrected age of , bridging the gap between kinematic ages ($11-19$ Myr) and other age-dating methods, such as isochrones and lithium depletion boundary ($20-26$ Myr). We explore several association size metrics that can track the spatial extent of  over time, and we determine that minimizing the variance along the heliocentric curvilinear coordinate $\xi^\prime$ (i.e., toward the Galactic Center) offers the least random and systematic errors, due to the wider $UVW$ space velocity dispersion of members of  along the $U$-axis, which tends to maximize the spatial growth of the association along the $\xi^\prime$-axis over time.

§ INTRODUCTION

Nearby young associations (NYAs) are sparse, coeval, gravitationally unbound stellar populations located in the solar neighborhood. Formed within a few million years of each other from the collapse of a single molecular cloud or cloud complex, members of these kinematic associations have similar Galactic positions and space velocities, and they share the same age and chemical composition due to their common formation history. However, their members span relatively large projected angular areas on the sky, as a result of the proximity and low spatial density of NYAs, which has made identification of NYAs challenging, especially prior to the  mission [Gaia Collaboration et al., 2016].

Such populations of nearby, age-calibrated stars are ideal laboratories to study the last stages of stellar and exoplanetary formation, and they provide strategic locations to search for isolated planetary-mass objects and for the direct imaging of giant exoplanets, thanks to the more favorable contrast between the host star and its companion, which remain relatively bright in the near-infrared (NIR) at such young ages. This contrast is even more advantageous for low-mass, late-type stars, which make up the majority of NYAs.

For up to a few hundred million years after their formation, members of an NYA retain similar kinematics, until gravitational perturbations cause their Galactic orbits to become randomized enough that they become indistinguishable from unrelated field stars. Using the full 6D kinematics of their members (i.e., their $XYZ$ Galactic positions and $UVW$ space velocities, where $U = dX/dt$, $V = dY/dt$, and $W = dZ/dt$), it is possible to compute backward Galactic orbits and trace members' trajectories back to the epoch when the NYA's spatial extent was minimal, which is assumed to coincide with the epoch of stellar formation. Thus, traceback analysis can provide a kinematic age estimate for members of an NYA, independent of stellar evolution models, unlike the more usual isochrone or lithium depletion boundary (LDB) methods.

One of the greatest difficulties in tracing back stellar trajectories is the need for precise astrometric and kinematic measurements to compute the full 6D kinematics of members of an NYA. The  Data Release 3 (DR3) data products provide an unmatched sample of precise parallax and proper motion measurements [Gaia Collaboration et al., 2022]. However, despite the increase of available measurements over the  Early Data Release 3 [Gaia Collaboration et al., 2021], absolute radial velocity measurements remain relatively inaccurate and represent, by far, the largest contribution to the total error in 6D kinematics during backward Galactic orbit integration, especially given how errors on Galactic position grow as stars are projected further back in time.

Jointly discovered by [Barrado y Navascués et al., 1999] and [Zuckerman et al., 2001], using data from the  catalog, the  () is one of the youngest and nearest known NYAs. Over 40 stars, located at an average distance $\sim 35$ pc, with a distance range of $9-72$ pc, with computed full 6D kinematics and visible signs of youth consistent with membership in , are considered to be bona fide members of the association, and more than a hundred additional candidate members are known, for which kinematic measurements for final membership determination are needed [Torres et al., 2006, Schlieder et al., 2010, Kiss et al., 2011, Schlieder et al., 2012, Malo et al., 2014, Riedel et al., 2014, Binks & Jeffries, 2016, Riedel et al., 2017, Gagné et al., 2018, Miret-Roig et al., 2020].

The youth and proximity of  make its members ideal targets for the search and characterization of exoplanets through direct imaging. For instance, its eponymous star, , has two known giant exoplanets, a debris disk, and exocomets [Lagrange et al., 2009, Chauvin et al., 2012, Kiefer et al., 2014, Lagrange et al., 2019]. PSO J318.5338–22.8603, a $6.5^{+1.3}_{-1.0}$  free-floating planetary-mass object of spectral type L7, was discovered using data from the Pan-STARRS survey [Liu et al., 2013]. With an $XYZ$ Galactic position and an $UVW$ space velocity compatible with membership in , its age is the same as the rest of the association.

However, traceback age estimates for  and other NYAs are inconsistent with other age-dating methods. Age estimates for  using the isochrones and LDB methods ($20-26$ Myr; Mentuch et al., 2008, Yee & Jensen, 2010, Malo et al., 2014, Mamajek & Bell, 2014, Binks & Jeffries, 2014, Bell et al., 2015, Galindo-Guil et al., 2022) are significantly older than most traceback age estimates ($11-13$ Myr; Ortega et al., 2002, Song et al., 2003, Ortega et al., 2004, Miret-Roig et al., 2018), although more recent kinematic approaches ($17-19$ Myr; Crundall et al., 2019, Miret-Roig et al., 2020) have addressed specific and distinct sources of bias, and have allowed the gap between isochrones and LDB methods to be bridged. [Crundall et al., 2019] used a forward modeling approach, which circumvents the bias on traceback ages due to measurement errors (see Section <ref>) that arises when computing backward Galactic orbits for individual members, whereas [Miret-Roig et al., 2020] used a more traditional traceback approach but minimized sample contamination through a rigorous selection process and used robust association size metrics to track the spatial extent of  over time.

In this work, we aim to correct for various sources of systematic errors in the computation of the traceback age of , and to determine whether this can further reconcile the tension between kinematic ages and other age-dating methods. Our approach uses a clean sample of bona fide members of , free of any kinematic outlier or unresolved multiple system (see Section <ref>), along with data from  data products and other radial velocity surveys, as the basis to compute backward Galactic orbits using their full 6D kinematics. We account for the bias on traceback age estimates due to measurement errors in astrometric and kinematic data (see Section <ref>) and two separate biases on radial velocity measurements, gravitational redshift and convective blueshift (see sections <ref> and <ref>), all of which tend to artificially push the epoch of minimal association size closer to the current-day epoch. We also test several association size metrics to evaluate the spatial extent of NYAs over time, in order to determine which provide the most accurate and reliable age estimates (see Section <ref>).

This study is structured as follows. First, in Section <ref>, we describe how we selected a clean, vetted sample of bona fide members of  and corrected for biases on radial velocity measurements. In Section <ref>, we describe the numerical method used to derive a kinematic age estimate of NYAs, and we apply it to simulated samples of stars in order to assess the precision and bias of every association size metric. In Section <ref>, we apply our method to our sample of members of , and we compare our results to previous age estimates for . Finally, we conclude our analysis in Section <ref>.

§ SAMPLE SELECTION

When assembling a sample with the aim to perform traceback analysis, candidates must be carefully vetted in order to minimize contamination by unrelated older stars that show a kinematic compatible with the NYA. Indeed, our simulations show that traceback ages are not only less precise with the addition of kinematic outliers, they are also biased toward younger ages because the positions of kinematic outliers do not converge at the epoch of stellar formation along with actual members of the NYA (see Section <ref>).

Therefore, in order to mitigate this issue, we first assembled a sample of candidate members of  from the literature [Zuckerman et al., 2001, Malo et al., 2013, Riedel et al., 2014, Miret-Roig et al., 2020] with available full 6D kinematics. Candidates were also vetted by the Bayesian analysis tool [ is documented as <http://www.exoplanetes.umontreal.ca/banyan/banyansigma.php>.] [Gagné et al., 2018], which utilizes astrometric and kinematic measurements to establish the membership probability of a star in 30 known NYAs, represented by 3D ellipsoidal models in $XYZ$ and $UVW$ space [Lee & Song, 2018, Gagné et al., 2018, Lee & Song, 2019]. In total,  stars, located at an average distance of $28.22 \pm 0.02$ pc and a distance range of $9.719-83.99$ pc, were identified as part of this sample, presented in Table <ref> and hereafter referred to as the input sample. Over half of these stars are low-mass M dwarfs, and most show signs of youth such as high X-ray, UV, or H$\alpha$ emission, fast rotation, or lithium emission lines. Stars within the input sample were further vetted in order to flag and exclude possible unresolved multiple systems (see Section <ref>) and kinematic outliers (see Section <ref>).

§.§ Stellar kinematics

We sourced most parallax, proper motion, and absolute radial velocity data from the  data products [Gaia Collaboration et al., 2022], while  data were used for very bright stars [Perryman et al., 1997], for which astrometry is sometimes more precise in  than . We also used other kinematic measurements compiled from the literature, and we used an error-weighted average to combine all available measurements. This allows us to obtain more precise and reliable measurements in order to better confirm the membership of members of  with  and compute more precise backward Galactic orbits for every member, which will in turn yield a more precise kinematic age estimate and help limit contamination by multiple systems and kinematic outliers. In addition, reliable and precise astrometric and kinematic measurements are also essential to minimize the bias on traceback ages due to measurement errors (see Section <ref>). The average $XYZ$ Galactic positions and $UVW$ space velocities of members of the input sample are presented in Table <ref>. The astrometric and kinematic measurements of all members of the input sample of  are presented in Table <ref>.

§.§ Unresolved multiple systems

A large fraction of stars in the Galaxy are in fact part of multiple systems, for which the added stellar motion around the barycenter of the system, along the line of sight, modulates radial velocity measurements. Unresolved multiple systems can therefore appear as false positives when one attempts to select NYA members based on kinematic data if a stars' $UVW$ space velocity happens to match the kinematic model of an NYA.

Fortunately, it is possible to flag these false positives by comparing radial velocity measurements at different epochs. Indeed, the radial velocity of stars in multiple systems changes periodically, and inconsistent measurements can thus be used as a way to flag unresolved multiple systems. Therefore, in order to exclude kinematic measurements contamination by the presence of these systems, we selected members from the input sample of  for which at least two radial velocity measurements are available in the literature and did not display variations $\geq 0.6$ , a value larger than the typical radial velocity scatter in our data set, indicating the possible presence of an unresolved companion. Also excluded were known spectral binaries from the literature and the Washington Double Star (WDS) catalog, stars with excess brightness in the color-magnitude diagram (CMD), as well as stars with a  Renormalized Unit Weight Error (RUWE) $\geq 2.0$.[The RUWE is documented in the Gaia Data Release Documentation at <https://gea.esac.esa.int/archive/documentation/GDR2/Gaia_archive/chap_datamodel/sec_dm_main_tables/ssec_dm_ruwe.html>.] This statistical indicator is a way to assess the reliability of astrometric data. It is expected to be $\sim 1.0$ for sources for which the single-star astrometric solution is a good fit [Lindegren, 2018] while a RUWE $> 1.4$ is a sign of an unresolved binary system [Stassun & Torres, 2021, El-Badry et al., 2021]. We chose a slightly less robust threshold because we also require stars to show relatively stable radial velocities over time. The various multiplicity indicators of members of the input sample of  are presented in Table <ref>.

The result of this selection process is a clean, uncontaminated subset of the input sample of  made up of  stars located at an average distance of $29.587 \pm 0.006$ pc and a distance range of $9.719-71.55$ pc, hereafter referred to as the core sample of , which can be used to compute the traceback age of the association. The average $XYZ$ Galactic positions and $UVW$ space velocities of members of the core sample are presented in Table <ref>. We note that all parameters of the input and core samples are similar, apart from the number of stars included in the sample and the $UVW$ space velocity dispersion ($\sigma_{UVW}$), which is significantly tighter for the core sample, as expected. The remaining  stars excluded from the core sample by this selection process, hereafter referred to as part of the extended sample, are not considered in the computation of the traceback age. The following stars deserve additional discussion:

HD 165189 This star is a slight kinematic outlier with only three available radial velocity measurements. It is more than 1 magnitude brighter than other members of  in the CMD, and the same is true of its companion located at a 1" separation. The absolute V vs. V–J is comparable to [Bell et al., 2015] and seems to best match the  parallax. Therefore, the  and EDR3 photometry might be unreliable. For now, it remains a   candidate and was kept in the extended sample.

HD 207043 This star's position in a  CMD is inconsistent with other members of , which puts into question its membership in the association. For now, it remains a  candidate and was kept in the extended sample.

AF Psc With updated kinematics, this triple star system with separations of 19" and 1000" [Malo et al., 2014, Kraus et al., 2014, Miret-Roig et al., 2020] is a poor kinematic match to the  model of . As a result, it was excluded from the extended sample.

G 271–110 With updated kinematics, this star is a poor kinematic match to the  model of , and its position in a CMD is also fainter than other members of . [Malo et al., 2014] noted that it also does not show a lithium absorption, whereas it normally should at its temperature if it were a member of . Also, at $1.71$, its  RUWE is a bit high, so it may be an unresolved binary system. As a result, it was excluded from the extended sample.

HD 14082 A This star displays slight radial velocity variations. It was excluded from the extended sample, and its co-mover, HD 14082 B, located at a separation of 14", was used instead.

HD 173167 This star is a spectroscopic binary system with few radial velocity measurements that display slight variations [Elliott & Bayo, 2016]. This star was excluded from the extended sample, and its co-mover, Smethells 20, was used instead.

§.§ Gravitational redshift

As a result of general relativity, radial velocity measurements are biased due to the gravitational redshift as light emitted at the surface of a star leaves its gravitational potential. Light is also later blueshifted slightly by the Sun's and the Earth's gravitational potentials, although this effect is negligible in comparison and safely ignored. The gravitational redshift is given as

\begin{equation}\label{eq:grav}
1 + z_{grav} = \frac{\lambda_{obs}}{\lambda} = {\left( 1 - \frac{2GM}{c^2 R} \right)}^{-\slfrac{1}{2}}.
\end{equation}

A star's radial velocity (RV) is obtained indirectly by measuring its Doppler shift, as it moves away from or toward the observer. It is estimated by its non-relativistic expression:

\begin{equation}\label{eq:doppler}
1 + z_{Doppler} = \frac{\lambda_{obs}}{\lambda} = 1 + \frac{RV}{c}.
\end{equation}

By combining equations <ref> and <ref>, the expected bias on traceback ages due to the gravitational redshift of a star can be estimated as

\begin{equation}\label{eq:rv_bias}
\Delta RV_{grav} = c \left( {\left( 1 - \frac{2GM}{c^2 R} \right)}^{-\slfrac{1}{2}} - 1 \right).
\end{equation}

Therefore, a star will appear to have a greater radial velocity (i.e., it will appear to move away faster from the observer). From equation <ref>, this bias can be derived by two fundamental parameters: the stellar mass ($M$) and radius ($R$). These values, however, cannot be measured directly for most stars, and an approximation based on spectral type is used instead (see Figure <ref>).
This approximation is calibrated on interferometric radius measurements of main-sequence stars [Ligi et al., 2016], semi-empirical spectral energy distribution (SED)-based determinations of angular radii of members of the core sample of  by [Pecaut & Mamajek, 2013], updated with  parallax measurements, as well as ensembles of dynamical masses of eclipsing binaries for young and pre-main-sequence stars compiled by [Hillenbrand & White, 2004] and other sources, updated with  parallax measurements when relevant. We avoid basing our approximation on absolute photometry or bolometric luminosity, due to the potentially large bias that we would obtain for yet unknown unresolved multiple systems in the core sample.

While the impact is relatively small,  for most stars, traceback analysis is highly sensitive to such biases: because members of the closest NYAs are spread across a large fraction of the sky, this bias will create the illusion that members are moving away from the observer in multiple directions, causing them to appear to reach a minimal spatial extent at a more recent epoch when performing a traceback analysis. More distant NYAs, which do not cover a large fraction of the sky, are proportionally less affected by this bias.

This bias also affects the determination of membership probabilities for NYAs with , because this process also utilizes radial velocities as one of its input parameters, although the impact is much less pronounced, given that the inherent velocity dispersion of young associations is larger, at $\sim 1-3$ .

§.§.§ Stellar masses

The mass of members of the core sample of  was estimated using their spectral type in order to avoid biases in cases of potential unresolved multiple systems that may remain in this sample. We used the [Pecaut & Mamajek, 2013] spectral type–mass relation[See <https://www.pas.rochester.edu/ emamajek/EEM_dwarf_UBVIJHK_colors_Teff.txt>.] (P13 relation hereafter) and compared it to dynamical mass measurements of known A2–M5 young stars ($2-150$ Myr) from the literature, improved with  parallax measurements when relevant. Figure <ref> shows that the dynamical masses of young stars are slightly larger than those of old field stars at a fixed spectral type. This is expected, given that stars warm up slightly in their pre-main-sequence phase, which lasts for $20-100$ Myr for $\leq 1$  stars (e.g., see Choi et al., 2016).

We built a young version of the P13 spectral type–mass relation by adjusting a two-segment linear fit to the difference in masses between the P13 relation and measured dynamical masses in log space to shift the P13 relation upward slightly, as shown in Figure <ref>. We found a relative standard deviation of $11\,\%$ for the dynamical masses of young stars with respect to our modified young sequence, which we adopt as $1 \sigma$ uncertainties on our spectral-type-based mass estimations.

§.§.§ Stellar radii

We similarly used the P13 spectral type–radius relations to determine the radii of the stars in the core sample based on their spectral types, to avoid potential biases due to unresolved binaries. We modified the main-sequence version of the P13 relations using the sample of known members of  in P13 to obtain a spectral type–radius relation that is appropriate for the younger age of , at a phase when stars are known to have an inflated radius due to their slow contraction during the pre-main-sequence phase (e.g., see Burrows et al., 1997), which is slowed down significantly in the case of low-mass stars by strong magnetic fields driven by their fast rotations (e.g., see Malo et al., 2014). P13 measured angular semi-diameters based on an SED fitting method [Masana et al., 2006], which we translated into stellar radii using the best available parallax measurements, mostly based on . We calculated a relative standard deviation of $9\,\%$ for the radius measurements with respect to our modified spectral type–radius sequence, which we adopt as our measurement error. We show the resulting sequence in Figure <ref>.

For both the masses and radii, we used the [Burrows et al., 1997] evolutionary models to extrapolate our empirical sequences beyond spectral type M5, and we assigned larger relative measurement errors of $20\,\%$ to account for possible model systematics, which may be especially important for young low-mass stars.

By combining the spectral type–mass and spectral type–radius relations described above, we built a spectral type–gravitational redshift relation shown in Figure <ref>. The application of this correction to radial velocity measurements of members of the core sample of  is presented in Table <ref>.

§.§ Convective blueshift

With the presence of convection cells at the stellar surface, light is emitted from areas moving toward or away from the observer, resulting in a Doppler broadening of spectral lines. Because rising gas is hotter and brighter, and it accounts for a larger fraction of the stellar surface, the broadening is uneven and causes a change in the shape of spectral lines and a net blueshift of the spectrum.

Convective blueshift has historically been difficult to measure. By comparing precise astrometric radial velocity measurements from the  and  missions with spectroscopic measurements, it is possible to isolate phenomena, other than the Doppler effect, that can shift spectral lines [Leão et al., 2019]. Convective blueshift can also be measured by studying the shape of the differential of spectral lines [Meunier et al., 2017]. Its impact on the instrumental point spread function is dependent on the instrumental resolution, the measurement method, and the spectral range. Lower stellar mass and higher stellar activity are linked to lower levels of convective blueshift, while lower metallicity and higher effective temperature are linked to higher levels. On the other hand, there is little effect from age or cyclic activity [Liebing et al., 2021]. The net radial velocity shifts derived by cross-correlation also depend on the spectral resolution of the instrument (i.e., the lower the resolution, the smaller the sensitivity to changes in the shapes of spectral lines) and the wavelength range [Allende Prieto et al., 2013].

The effect of convective blueshift is smaller than the aforementioned gravitational redshift, so the sum of both effects is expected to be a redshift (see Figure <ref>). For spectral types earlier than F2, we assume no convective blueshift, due to the absence of convection cells at the stellar surface. A three-segment polynomial sequence was fitted to blueshift measurements from the literature for F2–F7, F7–K4, and K4–M5 stars with a conservative estimated error of $0.2$ , to account for the dependence of the impact of convective blueshift on the instrumental point spread function. For stars later than M5, although they are expected to be fully convective, no robust measurements of convective blueshift are available yet. The case of YZ CMi [Baroch et al., 2020], a member of , is consistent with a null convective blueshift within the error bars. For this reason, we adopt a null correction with a $0.2$  measurement error until more data are available. The application of this correction to radial velocity measurements of members of the core sample of  is presented in Table <ref>.

§ TRACEBACK ANALYSIS

Traceback analysis is performed with the custom Python package [: Kinematic Age for Nearby Young Associations. The source code is available at <https://github.com/***/kanya>.], which is capable of computing the traceback age of a given sample of stars. The backward Galactic orbit of each star is computed and the spatial extent of the association is determined using a range of association size metrics (see Section <ref>) at every temporal step, in order to identify the epoch of minimal association size. As will be demonstrated in this section, this method is sensitive to sample contamination, biases introduced by measurement errors, and systematic effects that impact absolute radial velocity measurements. It is therefore crucial to address each of these potential issues in order to determine reliable traceback ages.

§.§ Backward orbital integration

A correction to radial velocity measurements was first applied before computing backward Galactic orbits, in order to compensate for the effects of gravitational redshift and convective blueshift (see sections <ref> and <ref>). Then the full current-day 6D kinematics of members of the core sample of  are computed (see Table <ref>). We use the galpy[galpy is documented at <http://github.com/jobovy/galpy>.] Python package [Bovy, 2015] with the Galactic potential model I from [Irrgang et al., 2013] to compute independent backward Galactic orbits for every member of the NYA. This model, which is an updated version of the Galactic potential from [Allen & Santillan, 1991], assumes a $9.5 \times 10^9$  spherical bulge, a $6.6 \times 10^{10}$  [Miyamoto & Nagai, 1975] disk, and a $1.8 \times 10^{12}$  spherical halo. This is the same Galactic potential used by [Miret-Roig et al., 2020], who concluded that variations in the kinematic age caused by the choice of the Galactic potential are smaller than the main source of uncertainty, due the short integration time for such young NYAs.

$XYZ$ Galactic coordinates are transformed from a heliocentric, right-handed Cartesian system into a Galactocentric, left-handed cylindrical system [Quillen et al., 2020]. The following values for the Sun's current-day Galactic position are adopted for this transformation:

\begin{align}\label{eq:sun_position}
R_\odot &= 8.12 \pm 0.03\textrm{\,kpc}\notag\\
\phi_\odot &= 0\textrm{\textdegree}\\
z_\odot &= 5.6 \pm 5.8\textrm{\,pc},\notag
\end{align}

where  is the Sun's Galactocentric radius [Gravity Collaboration et al., 2018] and $z_\odot$ is the Sun's height above or below the Galactic plane [Reid et al., 2019]. The Sun's Galactocentric longitude $\phi_\odot$ is null by definition. We adopt the following peculiar Solar motion to transform $UVW$ space velocities [Schönrich et al., 2010]:

\begin{align}\label{eq:sun_peculiar_velocity}
U_\odot &= 11.1_{-0.8}^{+0.7}\textrm{\,\kms}\notag\\
V_\odot &= 12.2 \pm 0.5\textrm{\,\kms}\\
W_\odot &= 7.3 \pm 0.4\textrm{\,\kms}.\notag
\end{align}

We also adopt the following local standard of rest (LSR) rotational velocity [Schönrich et al., 2010]:

\begin{equation}
V_{\textrm{LSR}} = 233 \pm 1.4 \textrm{\,\kms}.
\end{equation}

Once Galactic orbits are computed for every member of the NYA, their individual Galactocentric positions and velocities are transformed back into $XYZ$ Galactic coordinates and $UVW$ space velocities.

Though this approach is more physically accurate than simply ignoring the Galactic potential's impact on past trajectories, its inclusion may have little impact on both the bias and error of our method, for an association as young as $\sim 24$ Myr (see Figure <ref>). Members of NYAs have similar $XYZ$ Galactic positions and $UVW$ space velocities, and thus they follow similar Galactic orbits as well. Because only the members' relative positions to the average position of the NYA are considered, the Galactic potential may have little impact on traceback ages.

§.§ Kinematic outliers

Stars in the core sample (see Section <ref>) were further investigated in order to identify kinematic outliers. First, the $XYZ$ Galactic position and $UVW$ space velocity standard deviations are computed independently along each 6D component at every temporal step of the stars' trajectory. Stars that stray away from the core of the NYA, beyond a $3\sigma$ threshold in either $XYZ$ Galactic position or $UVW$ space velocity at any point and in any component, are flagged as kinematic outliers and excluded from the computation of size metrics. This process is repeated recursively until no more stars are flagged as kinematic outliers.

The robust covariance estimator developed by [Pedregosa et al., 2011] as part of the scikit-learn[Scikit-learn's robust covariance estimator is documented at <https://scikit-learn.org/stable/modules/covariance.html?highlight=robust%20covariance#robust-covariance-estimation>.] Python package was then used to identify kinematic outliers independently at every temporal step of the traceback. Stars identified as outliers at least $70\,\%$ of the time along their trajectory were flagged as kinematic outliers and ignored in the computation of all association size metrics. Scikit-learn's robust covariance estimator uses the empirical covariance matrix ($\Sigma$) of the stars' $XYZ$ Galactic positions and $UVW$ space velocities:

\begin{equation}
\label{eq:covariance_matrix}
\Sigma = E \left( (\boldsymbol{x} - E(\boldsymbol{x})) (\boldsymbol{x} - E(\boldsymbol{x}))^T \right),
\end{equation}

where $E(\boldsymbol{x})$ is the expected value of a random vector $\boldsymbol{x}$. Since the empirical covariance matrix is sensitive to the presence of outliers in the data, it is necessary to perform outlier detection in order to compute a robust estimator of the covariance matrix. This is accomplished by computing the Mahalanobis distance [Mahalanobis, 1936] of each star. For an observation $x_i$ of a distribution with an empirical covariance matrix $\Sigma$ and a mean $\mu$, the Mahalanobis distance ($d_M$) to the mode is given by

\begin{equation}
\label{eq:Mahalanobis}
d_M^2(x_i) = (x_i - \mu)^T \Sigma^{-1} (x_i - \mu).
\end{equation}

This approach is an effective tool to perform outlier detection in noisy and irregular data sets because it takes into account covariances in the distribution. In summary, the distance is measured and scaled to the standard deviation along each principal component of the 6D distribution.

The number of stars included in the sample is a limitation to the precision and accuracy of the traceback method. As one would expect, traceback ages are more precise and accurate with a larger sample. However, as more stars are included in the sample, some will likely be kinematic outliers or unrecognized binary stars, and these will in turn both negatively affect the precision of traceback ages and create an additional bias toward younger ages. The trajectories of these stars tend to deviate significantly from the average trajectory of the association, and even just a few such outliers can dramatically lower the traceback age by artificially increasing the dispersion in $XYZ$ Galactic positions and $UVW$ space velocities. Therefore, it is advantageous to limit traceback samples to smaller ensembles of stars, free of such contamination. However, including more stars in the sample will decrease the minimal theoretical error of ages determined with the traceback method, because the epoch of minimal association size will be more clearly defined, provided that no kinematic outliers affect the results negatively.

§.§ Association size metrics

The empirical covariance matrices of members of an NYA were calculated at each temporal step in both $XYZ$ or  Galactic coordinate systems. The latter, identical to the system used by [Miret-Roig et al., 2020], is a heliocentric curvilinear coordinate system that minimizes variations along every individual component during the backward orbital integration by moving the origin on a circular orbit around the Galactic Center at a frequency of $\omega_\odot = V_{\textrm{LSR}} / R_\odot = 28.7 \textrm{\,km} \textrm{\,s}^{-1} \textrm{\,kpc}^{-1}$. In addition to the individual diagonal terms of the empirical covariance matrices, the determinant and the trace of the associated covariance matrices were investigated as plausible association size metrics, similarly to the approach taken by [Miret-Roig et al., 2020].

We also considered the spatial–kinematic cross-covariance terms for members at each temporal step as potential association size metrics, in both $XYZ$ and  Galactic coordinate systems. Covariances between spatial and kinematic terms of a given direction are expected to be minimal at the epoch of stellar formation and naturally increase over time (e.g., see Crundall et al., 2019). Once again, the determinant, trace, and individual diagonal terms of the cross-covariance matrices were included in the size metrics under consideration.

We implemented both “empirical” and “robust” versions of each association size metric based on the empirical $XYZ$ and  covariance matrices (equation <ref>). While empirical metrics give each member an equal weight, robust metrics assign them a weight inversely proportional to their Mahalanobis distance (see equation <ref>) to the association's core.

We also investigated the use of median absolute deviations (MADs) as an alternate set of association size metrics, given that the MAD depends more weakly on a small number of outliers. We computed the MAD along each spatial component, in both $XYZ$ and  Galactic coordinate systems, as a way to represent the spatial extent of an NYA at each epoch. This metric is not only less sensitive to kinematic outliers than the standard deviation, it is also expected to be even less sensitive to outliers than “robust” metrics described above, because the latter still assign a small but nonzero weight to significant outliers.

Finally, we computed the minimum spanning tree (MST) of the position of members at each epoch in both $XYZ$ and  Galactic coordinate systems using a Kruskal algorithm [Kruskal, 1956]. MSTs are undirected graphs that link all vertices with the minimal possible total edge length, without allowing for loops. The mean edge length and the MAD of the edge lengths of the MST are used as association size metrics, which remain valid even in highly non-Gaussian geometries. This approach has the advantage of determining a characteristic size while being less sensitive to the shape of the NYA. A robust version of the mean edge length of the MST, using the same weights as other robust size metrics, was also explored.

§.§ Error on the traceback age

The traceback age is highly dependent on which stars are included in the computation of association size metrics. To assess this error on the traceback age due to the sensitivity to sample definition, we employed a jackknife Monte Carlo approach. We used $1000$ iterations, each made up of a $50\,\%$ fraction of randomly selected stars the sample, to measure the standard deviation on the traceback age for every size metric.

We also made use of a Monte Carlo approach to compute the error on the traceback age due to the measurement errors in the astrometric and kinematic data. Randomized Gaussian fluctuations with a standard deviation equal to the reported measurement errors were added to radial velocity, proper motion, and parallax measurements for $1000$ iterations. The total error on the traceback age was computed with a Monte Carlo approach including both the jackknife and measurement errors.

§.§ Simulated samples

We constructed simulated NYA samples, designed to best represent the core sample of current-day members  in order to test and characterize our method precision, accuracy, and sensitivity to several sources of bias.
A model star was initialized with an $XYZ$ Galactic position and a $UVW$ space velocity equal to the current-day average values of members of the core sample of  (see Table <ref>). This model star's trajectory was traced back $24$ Myr to the birth epoch, a length of time similar to recent age estimates of  using the isochrones and LDB methods [Malo et al., 2014]. Then, a randomized Gaussian sample of  synthetic stars, the same number as the number of stars in the core sample of , was created, with average $XYZ$ Galactic positions and $UVW$ space velocities equal to those of the model star at the birth epoch, an initial $XYZ$ Galactic position scatter of $3.0$ pc, and an initial $UVW$ space velocity scatter of $1.0$ , along all three axes. It is impossible to know for sure what the initial $XYZ$ Galactic position scatter of  was. However, the value used is similar to the current-day scatter along the $\zeta^\prime$-axis, an axis with very little growth over time. The initial $UVW$ space velocity scatter was approximated by the average current-day $UVW$ space velocity scatter of members of the core sample of , because we assume that $UVW$ space velocities are relatively unchanged since the epoch of stellar formation for such a young NYA.

We used the extreme deconvolution algorithm[The extreme deconvolution algorithm is documented at <https://github.com/jobovy/extreme-deconvolution>.] developed by [Bovy et al., 2011] to mitigate the impact of elongation of $UVW$ space velocities along the line of sight, due to larger errors on radial velocity measurements than proper motion measurements. This general algorithm can infer a $d$-dimensional distribution function from a set of heterogeneous and noisy samples, and it can treat uncertainties properly. When applied on the current-day $UVW$ space velocities of members of the core sample of  and their respective uncertainties (see Figure <ref>), the resulting ellipsoids offer a more accurate description of distribution of $UVW$ space velocities in the core sample of , which in turn is used to create more representative simulated samples.

The simulated sample were then projected $24$ Myr forward in time such that synthetic stars are located as near to the Sun as current-day members of . $XYZ$ Galactic positions and $UVW$ space velocities were transformed into observables, and Gaussian-like fluctuations were added to the true sky coordinates, parallax, radial velocities, proper motions of each synthetic star in order to simulate measurement errors. A radial velocity shift was also added to simulate the bias on traceback ages due to the gravitational redshift and convective blueshift (see sections <ref> and <ref>). The errors on astrometric and kinematic measurements, and the radial velocity shifts were drawn from the real set of individual members of the core sample of . The $24$ Myr old simulated, biased, and noisy sample is then traced back over $50$ Myr, and its size is computed with the full array of association size metrics (see Section <ref>) in order to find the epoch of minimal spatial extent.

§.§.§ Sensitivity to kinematic measurement errors

An important bias to account for when computing the age of a coeval group of stars with traceback methods is the impact of Gaussian-like measurement errors in astrometric and kinematic data that may not simply result in a Gaussian-like distribution in the estimated epoch of minimal association size, but rather in a more complex probability distribution that will be not only inflated by an additional error term but also systematically biased with respect to the true epoch of minimal size closer to the current-day epoch.

The total observed Galactic position scatter ($\sigma_{\textrm{observed}}$) measured by association size metrics (see Section <ref>) is made of two components: the inherent Galactic position scatter of members of an NYA ($\sigma_{\textrm{inherent}}$) and an additional, artificial scatter due to the impact of measurement errors ($\sigma_{\textrm{error}}$). Assuming Gaussian distributions, the total observed scatter is given by:

\begin{equation}
\label{eq:scatters}
\sigma^2_{\textrm{observed}}(t) = \sigma^2_{\textrm{inherent}}(t) + \sigma^2_{\textrm{error}}(t).
\end{equation}

Hence, the total observed scatter (the scatter that is directly observed in the traceback analysis) is an approximation of the inherent scatter (the scatter we are trying to measure). Unlike the inherent scatter, which we assume is minimal at the epoch of stellar formation and grows over time, the scatter due to measurement errors is minimal at the current-day epoch and increases as stellar trajectories are traced back in time, due to the increasingly imprecise $XYZ$ Galactic positions that result from errors on $UVW$ space velocities, mostly on account of relatively imprecise radial velocity measurements. Adding Gaussian-like fluctuations to the current-day radial velocities of members of the core sample of  across all directions on the sky affects the angle of convergence of the reconstructed $UVW$ space velocities, and thus biases the epoch and $XYZ$ Galactic position where members converge. Therefore, traceback ages become not only less precise but also biased toward younger ages. The older the NYA, the longer the trajectories of its members must be traced back in time to find the epoch of minimal association size, and the greater the error and bias on traceback ages, limiting the traceback method to younger associations.

In an effort to characterize this bias on traceback ages, we constructed simulated NYA samples initialized with a range of initial $XYZ$ Galactic position and $UVW$ space velocity scatters. After initializing a set of synthetic samples following these distributions, we projected them forward in time and added a range of Gaussian-like fluctuations on the simulated radial velocity measurements and a range of radial velocity shifts (see Section <ref>). Then, the resulting biased, noisy, reconstructed $XYZ$ Galactic positions and $UVW$ space velocities were traced back in time to determine how the epoch of minimum association size was affected. The biases (i.e., the average offsets of traceback ages from the actual age) and errors (i.e., the width of the distributions of traceback ages) on the traceback ages that result from this analysis are presented in Figure <ref>.

Larger measurement errors on radial velocity result in accordingly larger biases and errors on the traceback age, due to the added, artificial $XYZ$ Galactic position scatter. To account for this bias, a correction to the traceback age is applied based on results from the simulation-derived bias that we observe for each association size metric, given that different size metrics may be more or less sensitive to this bias. These corrections are approximate, and they require assumptions on the approximate age of the association, its initial $XYZ$ Galactic position scatter, and its initial $UVW$ space velocity scatter. However, the resulting bias does not depend strongly on these assumptions as long as the initial $UVW$ space velocity scatter is not significantly smaller than $\sim 0.5$ (see Figure <ref>), in which case the traceback method becomes much more sensitive to this bias. A greater $UVW$ space velocity scatter is beneficial to the precision and accuracy of traceback ages, as the point of origin becomes more clearly defined, and the divergence of stars over time is faster and more easily measurable. In contrast, the initial $XYZ$ Galactic position scatter has little impact on the bias and error on the traceback age, although a smaller initial $XYZ$ Galactic position scatter does slightly improve the error on the traceback age. This can be understood from the simple fact that the inherent $XYZ$ scatter of low-$UVW$-scatter associations would grow much more slowly over time, reducing the signal that we rely on to determine the traceback age.

In the extreme case of a null $UVW$ scatter, the NYA would not grow over time and the only change in the spatial extent of the NYA would be due to the Gaussian-like fluctuations added to the simulated current-day observations, which reduces the error on the traceback age because the bias on the traceback age is close to the actual age of the simulated NYA (i.e., the spatial extent is always minimal at the current-day epoch). In fact, with greater artificial errors on observations with a null $UVW$ space velocity scatter, the result is even more extreme: a bias on the traceback age equal to the actual age of the simulated NYA, and a null error on the traceback age. While this case is not physically realistic, it highlights how low-$UVW$-scatter NYAs are more sensitive to the bias on the traceback age due to measurement errors.

In order to compute a correction for the bias due to measurement errors in our analysis, we created simulated samples with the same $XYZ$ Galactic positions and $UVW$ space velocities, initial scatters, and measurement errors, the same number of members and the same age as those used in for simulated samples in Figure <ref>, and a null radial velocity shift, as this bias is corrected independently (see Section <ref>). With these assumptions, we computed a  correction to the computed traceback age of , depending on the size metric to account for the bias due to measurement errors.

§.§.§ Sensitivity to radial velocity shift

Another important source of bias is the shift on radial velocity measurements due to both gravitational redshift and convective blueshift, which add a spectral-type-dependent shift of  (see sections <ref> and <ref>).

Our simulated samples (see Figure <ref>) show that this can bias the computed traceback age of members of the core sample of , for which the typical range of total radial velocity shift is $0.2-0.6$ , by $\sim 4.0$ Myr. A greater bias on radial velocity increases the bias on the traceback age, although the error on the traceback age decreases slightly. Therefore, it is necessary to account for the effects of the gravitational redshift and convective blueshift on radial velocity measurements in order to derive reliable traceback ages. In order to compensate for this bias, a spectral-type-dependent correction equal to the opposite of the computed total radial velocity shift (see Figure <ref>), is applied to empirical radial velocity measurements of members of the core sample of .

§ RESULTS AND DISCUSSION

We applied the traceback method described in Section <ref> to the core sample of  members described in Section <ref> to derive a kinematic age estimate for the association.

§.§ Traceback analysis of the

The traceback analysis of , back to $50$ Myr into the past, is presented in figures <ref> and <ref>. With  data, one member included in the core sample, CD–31 16041, was automatically identified as a kinematic outlier (see Section <ref>), due to its Galactic position reaching $4.3\sigma$ above the association's average Galactic position along the $\xi^\prime$-axis, beyond the $3\sigma$ threshold (se Section <ref>). Its RUWE does not suggest it is unresolved multiple system. However, with a total of 12 available radial velocity measurements from  and other kinematic surveys, this star is no longer considered an outlier and it was used in the computation of association size metrics.

As expected from their similar kinematics, members of the core sample of  follow similar backward Galactic orbits. In Figure <ref>, the dispersion along the $\xi^\prime$-axis reaches a clearly defined minimum value, which corroborates our assumption that the wider $UVW$ space velocity dispersion of members of  along the $U$-axis would result a greater change in association size over time. In contrast, no clear minimum is observed along the $\eta^\prime$- and $\zeta^\prime$-axes. As expected, members follow sinusoidal trajectories along the $\zeta^\prime$-axis (perpendicular to the Galactic plane) and do not reach a minimum at the epoch of stellar formation. Because members do not clearly converge across these two directions, there is little useful data along these two axes for the computation of the traceback age of .

§.§ Traceback age of the

Figure <ref> shows several association size metrics back to $35$ Myr in the past for members of the core sample of . The traceback ages of  for all association size metrics tested in this work are presented in Table <ref>. These results confirm that, due to the wider $UVW$ space velocity dispersion of members of  along the $U$-axis ($\sigma_U = 1.66$ ), size metrics along the $X$- and $\xi^\prime$-axes show greater contrast (i.e., the relative change from the current-day epoch to the epoch of minimal association size), and reach a minimum value at a more distant epoch when compared to size metrics along the $Y$- ($\sigma_V = 0.52$ ), $\eta^\prime$-, $Z$- ($\sigma_W = 0.73$ ) or $\zeta^\prime$-axes, consistent with the result from our simulated samples (see Section <ref> and  Figure <ref>). Size metrics along the $Z$- and $\zeta^\prime$-axes are safely ignored, due to the lack of convergence along these axes. Compound size metrics, which include data along all axes, such as the determinant and trace of the  covariance matrix, or the total  MAD, reach a minimum value at a slightly older epoch but with an inferior contrast as a result of the narrower $UVW$ space velocity dispersion of members of the core sample of  along the $\eta^\prime$- and $\zeta^\prime$-axes. We also note that size metrics along the $\xi^\prime$-axis have a greater contrast than those based on the $X$-axis. As a result, we focus the remainder of our analysis on size metrics along the $\xi^\prime$-axis.

The contrast is slightly higher in the case of robust size metrics than equivalent empirical size metrics, but the epoch of minimal association size is similar and their jackknife error is significantly larger. This is an expected result for robust size metrics, due to the smaller number of members used by these metrics. A similar observation can be made for size metrics based on the spatial–kinematic cross-covariance matrices such as the $X-U$ and $\xi^\prime-\dot{\xi}^\prime$ cross-variances: a minimal value is reached at the same epoch as for their counterparts based on the $XYZ$ and  covariance matrices, and the contrast is the highest of any size metric tested in this study. However, the jackknife error of $\xi^\prime-\dot{\xi}^\prime$ cross-variance is much larger, suggesting that it is unsuitable for this traceback analysis.

Size metrics based on the MAD offer contrasts similar to those of their counterparts based on the spatial–kinematic cross-covariance matrices. However, their total errors are larger than those of size metrics based on the empirical $XYZ$ and  covariance matrices. As for size metrics based on the use of an MST, we observe no difference with respect to those based on the $XYZ$ or  coordinate systems. This result is not too surprising, given how MSTs are less sensitive to the shape of the NYA and only measure the change in distances between members. Nevertheless, size metrics based on the use of an MST are excluded from the current analysis, due to their low contrast.

We also tested the impact of taking into account the Galactic potential, correcting the radial velocity shifts, or minimizing the impact of the sample contamination by kinematic outliers and multiple systems on the computation of accurate traceback ages for a young NYA like . We find traceback ages  younger without considering the impact of the Galactic potential. As expected based on our discussion in Section <ref>, the impact of the Galactic potential on traceback ages is small, but not insignificant. The impact of the choice of the Galactic potential is negligible: if we use the Galactic potential from [Bovy, 2015] instead of model I from [Irrgang et al., 2013], the final age difference is $< 0.2$ Myr, confirming the result from [Miret-Roig et al., 2020]. If radial velocity measurements are not corrected for, the unaccounted impacts of gravitational redshift and convective blueshift cause our traceback ages to be  younger. The biggest change, however, comes from mitigating contamination by kinematic outliers and multiple systems. With the full input sample, traceback ages for  are  younger.

We found a raw traceback age of   for  using the $\xi^\prime$ variance as the size metric. We computed a correction to account for the bias due to measurement errors for every size metric with simulated samples representative of  using the same parameters as Figure <ref> to account for the limited number of members and the impact of measurement errors in astrometric and kinematic data (see Section <ref>). The measurement errors that were applied on simulated stars were taken from one-to-one associations between simulated stars and real members of the core sample of . We find that the measurement errors specific to our sample alone cause a  bias toward younger ages, using the $\xi^\prime$ variance as the size metric. Therefore, we apply a $0.6$ Myr correction to our raw traceback age and find a final, corrected traceback of  for .

§.§ Comparison with other age-dating methods

While our corrected traceback age estimate for  is significantly older than previous traceback age estimates ($11-13$ Myr; Ortega et al., 2002, Song et al., 2003, Ortega et al., 2004, Miret-Roig et al., 2018), it is compatible with several more recent studies. Figure <ref> shows the age distribution of  using the corrected $\xi^\prime$ variance as the association size metric, compared to other recent kinematic age estimates for  by [Crundall et al., 2019] and [Miret-Roig et al., 2020]. All three results are compatible with each other despite the differences in samples and methods.

The result from [Crundall et al., 2019], who reported a kinematic age of  for , can be more directly compared to our results using the $X-U$ spatial–kinematic cross-covariance matrix, because this is the same metric indirectly used by the  method (the age of their forward model is constrained by the assumption that the initial spatial–kinematic model of  has no spatial–kinematic cross-covariances). Using this size metric, we report a slightly older traceback age of $19.9 \pm 3.1$ Myr, including a $1.0$ Myr correction to account for the biases on traceback ages due to measurement errors. This difference is mainly due to the more reliable radial velocity measurements from  and other kinematic surveys used in our study.

[Miret-Roig et al., 2020] reported traceback ages of  and  for  by minimizing the trace and the determinant of the  covariance matrix, respectively, computed with scikit-learn's robust estimator [Pedregosa et al., 2011]. This is compatible with our result using the same association size metrics. We report corrected traceback ages of $18.6 \pm 4.0$ Myr and $21.3 \pm 6.5$ Myr using the same size metrics, respectively, although the errors on the traceback age are significantly larger. If instead we use the exact same sample of $26$ bona fide members as described in [Miret-Roig et al., 2020], along with our robust version of the trace and determinant of the  covariance matrix, we find corrected traceback ages of $21.3 \pm 8.4$ Myr and $22.4 \pm 8.1$ Myr, respectively. However, in our analysis, we excluded these specific size metrics in favor of size metrics with greater contrast and lower jackknife errors, such as the $\xi^\prime$ variance. The difference between our result and that of [Miret-Roig et al., 2020] is likely due not to differences in samples, but rather to the methods that were used to compute the spatial extent of the association and the Galactic orbits of its members.

Nevertheless, all of the ages measured by traceback analyses remain younger than the age estimates of  determined using other approaches, such as the isochrones and LDB methods ($20-26$ Myr; Mentuch et al., 2008, Yee & Jensen, 2010, Malo et al., 2014, Mamajek & Bell, 2014, Binks & Jeffries, 2014, Bell et al., 2015, Galindo-Guil et al., 2022). While precise astrometric and kinematic measurements, clean uncontaminated samples, and corrections applied to radial velocity measurements to account for the gravitational redshift and convective blueshift helped bridge the gap with other methods, kinematic ages still remain younger than those calculated using the isochrones and LDB methods by $\sim 1-3$ Myr.

This difference may be explained by several factors. First, there may still be contamination in the core sample of , or a more robust way of measuring the spatial extent of the association over time might be necessary. However, the issue could also be more fundamental. Based on ages measured for star-forming regions such as Tau-Aur [Krolikowski et al., 2021], members may remain gravitationally bound to interstellar dust and gas for less than $\sim 2-3$ Myr after their formation. This would bias traceback ages, regardless of the association size metric used, because the initial assumption that members follow independent Galactic orbits would only be true once the initial cloud has been dissipated and members are free of the influence of the rest of the association. In other words, the traceback method would measure the time since members of the NYA are gravitationally unbound, not the time since stellar formation. Such a systematic difference could potentially be further tested if the greater kinematic accuracy allows us to measure the traceback ages of older associations in the near future, as we may expect to find a similar offset between traceback ages and other methods.

§ CONCLUSIONS

In this study, we created a numerical tool capable of the deriving the age of an NYA, based solely on astrometric and kinematic measurements by tracing back the Galactic orbits of a sample of stars and evaluating the size of the association with multiple metrics. The  Python package takes into account several observational biases, including the bias on traceback ages due to measurement errors in astrometric and kinematic data, as well as the biases on absolute radial velocity measurements due to the gravitational redshift and convective blueshift of the star. Our results confirm that it is crucial to compensate for these biases and limit sample contamination by unresolved multiple systems and other kinematic outliers, in order to compute reliable traceback age estimates.

We applied our method to a core sample of  members of , which were assembled from the literature and include data from the  catalog, and we found that minimizing the variance along the $\xi^\prime$-axis offers the least random and systematic errors, due to the wider $UVW$ space velocity dispersion of members of  along the $U$-axis, which tends to maximize its spatial growth along the $\xi^\prime$-axis over time. We found a corrected traceback age of , a result compatible with other recent kinematic ages found by other studies, but still slightly lower, by $\sim 1-3$ Myr, than ages obtained using either the isochrones or LDB methods.

In future works, we plan to apply our traceback method to other known NYAs. Results will also be improved by new data from dedicated radial velocity surveys and future releases of the  catalog, which has dramatically increased the number of stars with precise radial velocity measurements. More robust NYA modeling and association size metrics may also help further bridge the gap with ages derived using non-kinematic methods.

§ ACKNOWLEDGMENTS

We would like to thank Miret-Roig, N. and Mamajek, E. for their assistance, useful discussions, and helpful answers to our inquires. J. G. and R. D. acknowledge the support of the Natural Sciences and Engineering Research Council of Canada (NSERC), funding reference numbers RGPIN-2021-03121 and RGPIN-2017-06777, respectively. This work was partially carried under a Banting grant from NSERC. This work has made use of: the SIMBAD database and VizieR catalog access tool, operated at the Centre de Données astronomiques de Strasbourg, France [Ochsenbein et al., 2000]; data products from the Two Micron All Sky Survey (2MASS), which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center (IPAC)/California Institute of Technology (Caltech), funded by the National Aeronautics and Space Administration (NASA), and the National Science Foundation (Skrutskie et al., 2006; DOI: 10.26131/IRSA2); data from the European Space Agency (ESA) mission
(Gaia Collaboration et al., 2016). Gaia data are being processed by the  Data Processing and Analysis Consortium (DPAC;
<https://www.cosmos.esa.int/web/gaia/dpac/consortium>). Funding for the DPAC
has been provided by national institutions, in particular the institutions
participating in the  Multilateral Agreement (MLA). The  mission website is <https://www.cosmos.esa.int/gaia>. The  archive website is <https://archives.esac.esa.int/gaia>.

D.C. wrote the codes and the manuscript, and generated figures and tables; J.G. compiled the list of candidates and built figures <ref>–<ref>; R.D. led the analysis and provided general comments.

[Gagné et al., 2018], galpy [Bovy, 2015], Scikit-learn [Pedregosa et al., 2011], Extreme deconvolution [Bovy et al., 2011], Astropy [Astropy Collaboration et al., 2013, Astropy Collaboration et al., 2018, Astropy Collaboration et al., 2022]

[Allen & Santillan, 1991]
Allen, C., & Santillan, A. 1991, , 22, 255

[Allende Prieto et al., 2013]
Allende Prieto, C., Koesterke, L., Ludwig, H. G., Freytag, B., &
Caffau, E. 2013, , 550, A103

[Allers et al., 2016]
Allers, K. N., Gallimore, J. F., Liu, M. C., & Dupuy, T. J. 2016, , 819,

[Allers & Liu, 2013]
Allers, K. N., & Liu, M. C. 2013, , 772, 79

[Astropy Collaboration et al., 2013]
Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013,
, 558, A33

[Astropy Collaboration et al., 2018]
Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al.
2018, , 156, 123

[Astropy Collaboration et al., 2022]
Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022,
, 935, 167

[Azulay et al., 2017]
Azulay, R., Guirado, J. C., Marcaide, J. M., et al. 2017, , 607,

[Baroch et al., 2020]
Baroch, D., Morales, J. C., Ribas, I., et al. 2020, , 641, A69

[Barrado y Navascués et al., 1999]
Barrado y Navascués, D., Stauffer, J. R., Song, I., & Caillault,
J. P. 1999, , 520, L123

[Bell et al., 2015]
Bell, C. P. M., Mamajek, E. E., & Naylor, T. 2015, , 454, 593

[Binks & Jeffries, 2014]
Binks, A. S., & Jeffries, R. D. 2014, MNRAS: Letters, 438, L11

[Binks & Jeffries, 2016]
Binks, A. S., & Jeffries, R. D. 2016, , 455, 3345

[Bovy, 2015]
Bovy, J. 2015, , 216, 29

[Bovy et al., 2011]
Bovy, J., Hogg, D. W., & Roweis, S. T. 2011, Annals of Applied
Statistics, 5, 1657

[Braun et al., 2021]
Braun, T. A. M., Yen, H.-W., Koch, P. M., et al. 2021, , 908, 46

[Bruntt et al., 2008]
Bruntt, H., North, J. R., Cunha, M., et al. 2008, , 386, 2039

[Burrows et al., 1997]
Burrows, A., Marley, M., Hubbard, W. B., et al. 1997, , 491, 856

[Chauvin et al., 2012]
Chauvin, G., Lagrange, A. M., Beust, H., et al. 2012, , 542, A41

[Choi et al., 2016]
Choi, J., Dotter, A., Conroy, C., et al. 2016, , 823, 102

[Corbally, 1984]
Corbally, C. J. 1984, , 55, 657

[Crundall et al., 2019]
Crundall, T. D., Ireland, M. J., Krumholz, M. R., et al. 2019, ,
489, 3625

[da Silva et al., 2009]
da Silva, L., Torres, C. A. O., de La Reza, R., et al. 2009, , 508,

[Dai et al., 2019]
Dai, D.-C., Li, Z., & Stojkovic, D. 2019, , 871, 119

[Desidera et al., 2015]
Desidera, S., Covino, E., Messina, S., et al. 2015, , 573, A126

[Durkan et al., 2018]
Durkan, S., Janson, M., Ciceri, S., et al. 2018, , 618, A5

[Egret et al., 1992]
Egret, D., Didelon, P., McLean, B. J., Russell, J. L., & Turon, C.
1992, , 258, 217

[El-Badry et al., 2021]
El-Badry, K., Rix, H.-W., & Heintz, T. M. 2021, , 506, 2269

[Elliott & Bayo, 2016]
Elliott, P., & Bayo, A. 2016, , 459, 4499

[Elliott et al., 2014]
Elliott, P., Bayo, A., Melo, C. H. F., et al. 2014, , 568, A26

[Faherty et al., 2016]
Faherty, J. K., Riedel, A. R., Cruz, K. L., et al. 2016, , 225, 10

[Gagné & Faherty, 2018]
Gagné, J., & Faherty, J. K. 2018, , 862, 138

[Gagné et al., 2014]
Gagné, J., Lafrenière, D., Doyon, R., Malo, L., & Artigau,
É. 2014, , 783, 121

[Gagné et al., 2018]
Gagné, J., Roy-Loubier, O., Faherty, J. K., Doyon, R., & Malo,
L. 2018a, , 860, 43

[Gagné et al., 2018]
Gagné, J., Mamajek, E. E., Malo, L., et al. 2018b,
, 856, 23

[Gaia Collaboration et al., 2016]
Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016,
, 595, A1

[Gaia Collaboration et al., 2018]
Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, ,
616, A1

[Gaia Collaboration et al., 2021]
—. 2021, , 649, A1

[Gaia Collaboration et al., 2022]
Gaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2022, arXiv
e-prints, arXiv:2208.00211

[Galindo-Guil et al., 2022]
Galindo-Guil, F. J., Barrado, D., Bouy, H., et al. 2022, , 664, A70

[Gizis et al., 2002]
Gizis, J. E., Reid, I. N., & Hawley, S. L. 2002, , 123, 3356

[Gravity Collaboration et al., 2018]
Gravity Collaboration, Abuter, R., Amorim, A., et al. 2018, , 615,

[Gray et al., 2006]
Gray, R. O., Corbally, C. J., Garrison, R. F., et al. 2006, , 132,

[Gray et al., 2003]
Gray, R. O., Corbally, C. J., Garrison, R. F., McFadden, M. T., &
Robinson, P. E. 2003, , 126, 2048

[Gunn et al., 1988]
Gunn, J. E., Griffin, R. F., Griffin, R. E. M., & Zimmerman, B. A.
1988, , 96, 198

[Hillenbrand & White, 2004]
Hillenbrand, L. A., & White, R. J. 2004, , 604, 741

[Houk, 1978]
Houk, N. 1978, Michigan Catalogue of Two-Dimensional Spectral Types for the
HD Stars (Ann Arbor: Dept. of Astronomy, University of Michigan)

[Houk & Smith-Moore, 1988]
Houk, N., & Smith-Moore, M. 1988, Michigan Catalogue of Two-dimensional
Spectral Types for the HD Stars. Volume 4, Declinations -26°.0 to
-12°.0., Vol. 4 (Ann Arbor: Dept. of Astronomy, University of

[Irrgang et al., 2013]
Irrgang, A., Wilcox, B., Tucker, E., & Schiefelbein, L. 2013, ,
549, A137

[Janson et al., 2018]
Janson, M., Durkan, S., Bonnefoy, M., et al. 2018, , 620, A33

[Jeffers et al., 2018]
Jeffers, S. V., Schöfer, P., Lamert, A., et al. 2018, , 614,

[Jönsson et al., 2020]
Jönsson, H., Holtzman, J. A., Allende Prieto, C., et al. 2020, ,
160, 120

[Katoh et al., 2013]
Katoh, N., Itoh, Y., Toyota, E., & Sato, B. 2013, , 145, 41

[Kervella et al., 2004]
Kervella, P., Thévenin, F., Morel, P., et al. 2004, in Stars as
Suns : Activity, Evolution and Planets, ed. A. K. Dupree & A. O. Benz,
Vol. 219, 80

[Kharchenko et al., 2007]
Kharchenko, N. V., Scholz, R. D., Piskunov, A. E., Röser, S., & Schilbach,
E. 2007, Astronomische Nachrichten, 328, 889

[Kiefer et al., 2014]
Kiefer, F., Lecavelier des Etangs, A., Boissier, J., et al. 2014, ,
514, 462

[Kiss et al., 2011]
Kiss, L. L., Moór, A., Szalai, T., et al. 2011, , 411, 117

[Kraus et al., 2014]
Kraus, A. L., Shkolnik, E. L., Allers, K. N., & Liu, M. C. 2014, ,
147, 146

[Krolikowski et al., 2021]
Krolikowski, D. M., Kraus, A. L., & Rizzuto, A. C. 2021, , 162, 110

[Kruskal, 1956]
Kruskal, J. B. 1956, Proceedings of the American Mathematical Society, 7, 48

[Lafarga et al., 2020]
Lafarga, M., Ribas, I., Lovis, C., et al. 2020, , 636, A36

[Lagrange et al., 2009]
Lagrange, A. M., Gratadour, D., Chauvin, G., et al. 2009, , 493,

[Lagrange et al., 2019]
Lagrange, A. M., Meunier, N., Rubini, P., et al. 2019, NatAs, 3, 1135

[Leão et al., 2019]
Leão, I. C., Pasquini, L., Ludwig, H. G., & de Medeiros, J. R.
2019, , 483, 5026

[Lee & Song, 2018]
Lee, J., & Song, I. 2018, , 475, 2955

[Lee & Song, 2019]
—. 2019, , 486, 3434

[Liebing et al., 2021]
Liebing, F., Jeffers, S. V., Reiners, A., & Zechmeister, M. 2021,
, 654, A168

[Ligi et al., 2016]
Ligi, R., Creevey, O., Mourard, D., et al. 2016, , 586, A94

[Lindegren, 2018]
Lindegren, L. 2018, gAIA-C3-TN-LU-LL-124

[Liu et al., 2016]
Liu, M. C., Dupuy, T. J., & Allers, K. N. 2016, , 833, 96

[Liu et al., 2013]
Liu, M. C., Magnier, E. A., Deacon, N. R., et al. 2013, , 777, L20

[Loehner-Boettcher et al., 2019]
Loehner-Boettcher, J., Schmidt, W., Schlichenmaier, R., Steinmetz, T.,
& Holzwarth, R. 2019, yCat, J/A+A/624/A57

[Madsen & Gaensler, 2013]
Madsen, G. J., & Gaensler, B. M. 2013, , 209, 33

[Mahalanobis, 1936]
Mahalanobis, P. C. 1936, in Proceedings of the National Institute of Sciences
of India, Vol. 2, 49–55

[Malo et al., 2014]
Malo, L., Artigau, É., Doyon, R., et al. 2014a, ,
788, 81

[Malo et al., 2014]
Malo, L., Doyon, R., Feiden, G. A., et al. 2014b, ,
792, 37

[Malo et al., 2013]
Malo, L., Doyon, R., Lafrenière, D., et al. 2013, , 762, 88

[Mamajek & Bell, 2014]
Mamajek, E. E., & Bell, C. P. M. 2014, , 445, 2169

[Masana et al., 2006]
Masana, E., Jordi, C., & Ribas, I. 2006, , 450, 735

[Mason et al., 2022]
Mason, B. D., Wycoff, G. L., Hartkopf, W. I., Douglass, G. G., &
Worley, C. E. 2022, yCat, B/wds

[Mentuch et al., 2008]
Mentuch, E., Brandeker, A., van Kerkwijk, M. H., Jayawardhana, R., &
Hauschildt, P. H. 2008, , 689, 1127

[Messina et al., 2010]
Messina, S., Desidera, S., Turatto, M., Lanzafame, A. C., & Guinan,
E. F. 2010, , 520, A15

[Meunier et al., 2017]
Meunier, N., Mignon, L., & Lagrange, A. M. 2017, , 607, A124

[Miret-Roig et al., 2018]
Miret-Roig, N., Antoja, T., Romero-Gómez, M., & Figueras, F. 2018,
, 615, A51

[Miret-Roig et al., 2020]
Miret-Roig, N., Galli, P. A. B., Brandner, W., et al. 2020, , 642, A179

[Miyamoto & Nagai, 1975]
Miyamoto, M., & Nagai, R. 1975, , 27, 533

[Montes et al., 2001]
Montes, D., López-Santiago, J., Fernández-Figueroa, M. J., &
Gálvez, M. C. 2001, , 379, 976

[Montes et al., 2001]
Montes, D., López-Santiago, J., Gálvez, M., et al. 2001, , 328,

[Montet et al., 2015]
Montet, B. T., Bowler, B. P., Shkolnik, E. L., et al. 2015, , 813,

[Moór et al., 2013]
Moór, A., Szabó, G. M., Kiss, L. L., et al. 2013, , 435,

[Newton et al., 2014]
Newton, E. R., Charbonneau, D., Irwin, J., et al. 2014, , 147, 20

[Nielsen et al., 2016]
Nielsen, E. L., De Rosa, R. J., Wang, J., et al. 2016, , 152, 175

[Nordström et al., 2004]
Nordström, B., Mayor, M., Andersen, J., et al. 2004, , 418, 989

[Ochsenbein et al., 2000]
Ochsenbein, F., Bauer, P., & Marcout, J. 2000, , 143, 23

[Ortega et al., 2002]
Ortega, V. G., de La Reza, R., Jilinski, E., & Bazzanella, B. 2002, , 575,

[Ortega et al., 2004]
Ortega, V. G., de la Reza, R., Jilinski, E., & Bazzanella, B. 2004,
, 609, 243

[Pecaut & Mamajek, 2013]
Pecaut, M. J., & Mamajek, E. E. 2013, , 208, 9

[Pedregosa et al., 2011]
Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, Journal of
Machine Learning Research, 12, 2825

[Pegues et al., 2021]
Pegues, J., Czekala, I., Andrews, S. M., et al. 2021, , 908, 42

[Perryman et al., 1997]
Perryman, M. A. C., Lindegren, L., Kovalevsky, J., et al. 1997, ,
323, L49

[Quillen et al., 2020]
Quillen, A. C., Pettitt, A. R., Chakrabarti, S., et al. 2020, ,
499, 5623

[Reid et al., 2019]
Reid, M. J., Menten, K. M., Brunthaler, A., et al. 2019, , 885, 131

[Riedel et al., 2017]
Riedel, A. R., Alam, M. K., Rice, E. L., Cruz, K. L., & Henry, T. J.
2017a, , 840, 87

[Riedel et al., 2017]
Riedel, A. R., Blunt, S. C., Lambrides, E. L., et al.
2017b, , 153, 95

[Riedel et al., 2018]
Riedel, A. R., Silverstein, M. L., Henry, T. J., et al. 2018, , 156,

[Riedel et al., 2014]
Riedel, A. R., Rodriguez, D., Finch, C. T., et al. 2014, The Astronomical
Journal, 147, 85

[Robertson & Hamilton, 1987]
Robertson, T. H., & Hamilton, J. E. 1987, , 93, 959

[Rodet et al., 2018]
Rodet, L., Bonnefoy, M., Durkan, S., et al. 2018, , 618, A23

[Schlieder et al., 2012]
Schlieder, J. E., Lépine, S., Rice, E., et al. 2012a,
, 143, 114

[Schlieder et al., 2010]
Schlieder, J. E., Lépine, S., & Simon, M. 2010, , 140, 119

[Schlieder et al., 2012]
—. 2012b, , 143, 80

[Schlieder et al., 2012]
—. 2012c, , 144, 109

[Schneider et al., 2019]
Schneider, A. C., Shkolnik, E. L., Allers, K. N., et al. 2019, ,
157, 234

[Schönrich et al., 2010]
Schönrich, R., Binney, J., & Dehnen, W. 2010, , 403, 1829

[Shkolnik et al., 2017]
Shkolnik, E. L., Allers, K. N., Kraus, A. L., Liu, M. C., & Flagg, L. 2017,
, 154, 69

[Shkolnik et al., 2012]
Shkolnik, E. L., Anglada-Escude, G., Liu, M. C., et al. 2012, , 758, 56

[Shkolnik et al., 2011]
Shkolnik, E. L., Liu, M. C., Reid, I. N., Dupuy, T., & Weinberger,
A. J. 2011, , 727, 6

[Simon & Schaefer, 2011]
Simon, M., & Schaefer, G. H. 2011, , 743, 158

[Simon et al., 2019]
Simon, M., Guilloteau, S., Beck, T. L., et al. 2019, , 884, 42

[Skrutskie et al., 2006]
Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, , 131,

[Song et al., 2003]
Song, I., Zuckerman, B., & Bessell, M. S. 2003, , 599, 342

[Song et al., 2012]
Song, I., Zuckerman, B., & Bessell, M. S. 2012, , 144, 8

[Soubiran et al., 2018]
Soubiran, C., Jasniewicz, G., Chemin, L., et al. 2018, , 616, A7

[Sperauskas et al., 2016]
Sperauskas, J., Bartašiūtė, S., Boyle, R. P., et al.
2016, , 596, A116

[Stassun & Torres, 2021]
Stassun, K. G., & Torres, G. 2021, , 907, L33

[Stephenson, 1986]
Stephenson, C. B. 1986, , 92, 139

[Teixeira et al., 2009]
Teixeira, R., Ducourant, C., Chauvin, G., et al. 2009, , 503, 281

[Torres et al., 2006]
Torres, C. A. O., Quast, G. R., da Silva, L., et al. 2006, , 460,

[Trifonov et al., 2020]
Trifonov, T., Tal-Or, L., Zechmeister, M., et al. 2020, , 636, A74

[Valenti & Fischer, 2005]
Valenti, J. A., & Fischer, D. A. 2005, , 159, 141

[van Leeuwen, 2007]
van Leeuwen, F. 2007, , 474, 653

[Weis, 1993]
Weis, E. W. 1993, , 105, 1962

[White et al., 2007]
White, R. J., Gabor, J. M., & Hillenbrand, L. A. 2007, , 133, 2524

[Yee & Jensen, 2010]
Yee, J. C., & Jensen, E. L. N. 2010, , 711, 303

[Zúñiga-Fernández et al., 2021]
Zúñiga-Fernández, S., Bayo, A., Elliott, P., et al. 2021,
, 645, A30

[Zuckerman et al., 2001]
Zuckerman, B., Song, I., Bessell, M. S., & Webb, R. A. 2001, , 562, L87