<figure><img src="image_1.png" /></figure>

# BANYAN. XI. The BANYAN Σ Multivariate Bayesian Algorithm to Identify Members of Young Associations with 150 pc

Jonathan Gagne<fnref n="1" />,<fnref n="8" /> ID, Eric E. Mamajek<fnref n="2" />,<fnref n="3" /> ID, Lison Malo<fnref n="4" /> ID, Adric Riedel<fnref n="5" /> ID, David Rodriguez<fnref n="5" /> ID, David Lafrenière<fnref n="4" /> ID, Jacqueline K. Faherty<fnref n="6" /> ID, Olivier Roy-Loubier<fnref n="4" /> ID, Laurent Pueyo<fnref n="5" /> ID, Annie C. Robin<fnref n="7" /> ID, and René Doyon<fnref n="4" /> ID

Received 2017 October 24; revised 2018 January 12; accepted 2018 January 26; published 2018 March 21

## Abstract

BANYAN Σ is a new Bayesian algorithm to identify members of young stellar associations within 150 pc of the Sun. It includes 27 young associations with ages in the range ~1–800 Myr, modeled with multivariate Gaussians in six-dimensional (6D) XYZUVW space. It is the first such multi-association classification tool to include the nearest sub-groups of the Sco-Cen OB star-forming region, the IC 2602, IC 2391, Pleiades and Platais 8 clusters, and the ρ Ophiuchi, Corona Australis, and Taurus star formation regions. A model of field stars is built from a mixture of multivariate Gaussians based on the Besançon Galactic model. The algorithm can derive membership probabilities for objects with only sky coordinates and proper motion, but can also include parallax and radial velocity measurements, as well as spectrophotometric distance constraints from sequences in color–magnitude or spectral type–magnitude diagrams. BANYAN Σ benefits from an analytical solution to the Bayesian marginalization integrals over unknown radial velocities and distances that makes it more accurate and significantly faster than its predecessor BANYAN II. A contamination versus hit rate analysis is presented and demonstrates that BANYAN Σ achieves a better classification performance than other moving group tools available in the literature, especially in terms of cross-contamination between young associations. An updated list of bona fide members in the 27 young associations, augmented by the Gaia-DR1 release, as well as all parameters for the 6D multivariate Gaussian models for each association and the Galactic field neighborhood within 300 pc are presented. This new tool will make it possible to analyze large data sets such as the upcoming Gaia-DR2 to identify new young stars. IDL and Python versions of BANYAN Σ are made available with this publication, and a more limited online web tool is available at http://www.exoplanetes.umontreal.ca/banyan/banyansigma.php.

**Key words:** brown dwarfs – methods: data analysis – proper motions – stars: kinematics and dynamics – stars: low-mass

**Supporting material:** figure set, machine-readable tables

## 1. Introduction

Coeval associations of stars that formed from a single molecular cloud are valuable benchmarks to study how the properties of stars evolve with time (e.g., Zuckerman & Song 2004; Torres et al. 2008). While precisely measuring the age of an individual star is challenging, the simultaneous study of a large ensemble of stars can provide age measurements with precisions down to a few Myr (e.g., Bell et al. 2015). A small number of young associations in the solar neighborhood that were identified to date are of particular interest in part because their lower-mass members can be studied easily. They are also coeval populations, making them valuable for measuring the initial mass function and serving as age calibrators for members across all masses.

Recent surveys (e.g., Gagné et al. 2014, 2015c, 2015b; Kellogg et al. 2015; Aller et al. 2016; Faherty et al. 2016; Liu et al. 2016) have started uncovering the substellar and planetary-mass members of nearby young stellar associations, which will make it possible to understand how their fundamental properties evolve with time.

Because young associations of the solar neighborhood are sparse and span up to ~20 pc in size, the distribution of their members can cover wide areas on the celestial sphere, and in some cases cover it almost entirely (e.g., the AB Doradus and β Pictoris moving groups; Zuckerman et al. 2001a, 2004). As a consequence, selection criteria based on sky coordinates and photometry alone are problematic.

As they formed recently and have not yet been perturbed significantly by other stars in the Galaxy, the members of a given young association still share similar space velocities UVW, with typical velocity dispersions below ~3 km s⁻<fnref n="1" />. This provides a way to identify the members of a young association; however, measuring their full kinematics requires not only proper motions, but also absolute radial velocities and parallaxes for every star. This constitutes the most challenging aspect in identifying their faint, low-mass members, as obtaining such measurements for a large sample of faint objects is prohibitive.

Various methods have been developed to identify members of young associations when only sky coordinates and proper motions are available. These include the convergent point tool (e.g., Mamajek 2005; Torres et al. 2006; Rodriguez et al. 2011), various goodness-of-fit metrics (e.g., Kraus et al. 2014;

<footnote n="1">¹ Carnegie Institution of Washington DTM, 5241 Broad Branch Road NW, Washington, DC 20015, USA; jgagne@carnegiescience.edu</footnote>

<footnote n="2">² Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA</footnote>

<footnote n="3">³ Department of Physics and Astronomy, University of Rochester, Rochester, NY 14627, USA</footnote>

<footnote n="4">⁴ 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</footnote>

<footnote n="5">⁵ Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA</footnote>

<footnote n="6">⁶ Department of Astrophysics, American Museum of Natural History, Central Park West at 79th St., New York, NY 10024, USA</footnote>

<footnote n="7">⁷ Institut UTINAM, OSU THETA, Univ. Bourgogne Franche-Comté, Besançon, France</footnote>

<footnote n="8">⁸ NASA Sagan Fellow.</footnote>
Bowler et al. 2017; Shkolnik et al. 2017), and the “good box” method (Zuckerman & Song 2004). Malo et al. (2013) developed BANYAN (Bayesian Analysis for Nearby Young AssociationNs), an algorithm based on Bayesian inference where moving groups are modeled with unidimensional Gaussian distributions in Galactic coordinates XYZ and space velocities UVW. More complex algorithms, such as BANYAN II (Gagné et al. 2014) and LACEwING (Riedel et al. 2017b) have more recently been developed, where associations are modeled with freely rotating tridimensional Gaussian ellipsoids in position and velocity spaces. BANYAN I included distance constraints from field and young sequences in a $M_J$ versus $I_C - J$ color-magnitude diagram in the 2MASS (Skrutskie et al. 2006) and Cousins (see Malo et al. 2013 for details) systems, and were defined for spectral classes K and M. BANYAN II included similar constraints based on two color-magnitude diagrams ($J - K_S$ versus $M_{W1}$ and $H - W2$ versus $M_{W1}$) in the 2MASS and WISE (Wright et al. 2010) systems, and were defined for spectral types later than M5.

These tools made it possible to identify hundreds of candidate members in nearby associations of stars, spanning the planetary to stellar-mass domains (e.g., Malo et al. 2013, 2014; Gagné et al. 2015b; Faherty et al. 2016; Liu et al. 2016). The majority of these classification tools include only the seven youngest (~10–200 Myr) and nearest (≤100 pc) moving groups, with the exception of LACEwING, which includes 16 associations and open clusters with a larger age range (~5–800 Myr).

The upcoming data releases of the Gaia mission (Gaia Collaboration et al. 2016b) will mark a new era in the study of young associations, as they will provide precise parallax measurements for a billion stars in the Galaxy, covering the full members of all associations within 150 pc down to late-M spectral types (Smart et al. 2017). This advancement in the census of members will improve kinematic models, lead to detailed measurements of initial mass functions, and open the door to the discovery of new sparse associations (e.g., see Oh et al. 2017). The first release of the Gaia mission (Gaia-DR1; Gaia Collaboration et al. 2016a) has already provided two million parallax measurements for stars in the Tycho-2 catalog (Høg et al. 2000), which have not yet been used to improve the membership classification tools described above.

This work presents BANYAN $\Sigma$, the next generation of the BANYAN tool based on Bayesian inference, which includes 27 associations with ages in the range ~1–800 Myr, completing the sample of known and well-defined associations within ~150 pc.<fnref n="9" /> BANYAN $\Sigma$ includes a significantly improved Gaussian mixture model of the Galactic disk which captures a larger fraction of field interlopers, and updated models of young associations that benefit from the most recent Gaia-DR1 parallax measurements. The models are also advanced to six-dimensional (6D) multivariate Gaussians that capture full correlations in the XYZUVW distribution of members, including those in mixed spatial-kinematic coordinates. Two versions of the BANYAN $\Sigma$ code (IDL<fnref n="10" /> and python) are made publicly available,<fnref n="11" /> and a web portal is made available for single-object queries.<fnref n="12" />

Most previously available classification tools rely on time-consuming algorithms, such as numerical integrals, which make it challenging to analyze large data sets such as the upcoming full Gaia release, and use various approximations in converting observables and kinematic models to probabilities that affect their classification performance. In BANYAN $\Sigma$, most of these approximations are removed, and Bayesian marginalization integrals are solved analytically. As a consequence, the tool is ≈80,000 times faster than its predecessor BANYAN II, making it easier to analyze very large data sets. BANYAN $\Sigma$ is also the first classification tool to include the Taurus, $\rho$ Ophiuchi and Corona Australis star-forming regions (e.g., Wichmann et al. 2000; Reipurth 2008), and the nearest OB association Sco-Cen, composed of the three subgroups: Upper Scorpius, Lower-Centaurus Crux, and Upper Centaurus Lupus (e.g., Blaauw 1946; de Zeeuw et al. 1999; Pecaut & Mamajek 2016). It is also the first such tool to include the IC 2602 (e.g., Whiteoak 1961; Mermilliod et al. 2009), IC 2391 (e.g., Platais et al. 2007), and Platais 8 (Platais et al. 1998) clusters, and one of the first to include the Pleiades cluster. For the last of these, Sarro et al. (2014) presented a multivariate Gaussian mixture model to assign membership probabilities based on kinematic and photometric observables. Their model uses a larger number of free parameters, made possible by the large number of known Pleiades members, but it does not include other young associations.

In Section 2, the framework of BANYAN II is described, which serves as a starting point for BANYAN $\Sigma$, described in detail in Section 3. An updated list of bona fide members for 27 young associations within 150 pc is presented in Section 4, and is used to build the multivariate Gaussian kinematic models of BANYAN $\Sigma$ in Section 5. Section 6 presents a multivariate Gaussian mixture model of field stars based on the Besançon Galactic model. A choice of Bayesian priors that ensures fixed recovery rates in all associations when using a $P = 90\%$ Bayesian probability threshold is described in Section 7. A performance analysis of BANYAN $\Sigma$ is presented in Section 8, and is compared to other tools available in the literature. The membership of stars previously considered as ambiguous are revisited in Section 9 using BANYAN $\Sigma$. This work is concluded in Section 10.

## 2. The Banyan II Algorithm

In this section, the Bayesian framework behind the BANYAN II tool (Gagné et al. 2014) is described. The framework of BANYAN $\Sigma$ will start from the same principles, but will include several improvements that are described in Section 3.

BANYAN II is a Bayesian classification algorithm, which uses the direct kinematic observables of a star {$O_i$}, namely its sky position ($\alpha, \delta$), proper motion ($\mu_\alpha, \mu_\delta$), radial velocity ($\nu$) and distance ($\varpi$), to determine the probability that the star belongs to a population described by a hypothesis $H_k$, corresponding to either the Galactic field or one of several young associations. This is done by applying Bayes' theorem:

$$P(H_k|\{O_i\}) = \frac{P(H_k)P(\{O_i\}|H_k)}{P(\{O_i\})}, \quad (1)$$

where the likelihood $P(\{O_i\}|H_k)$ is the probability that a member of $H_k$ displays the observables {$O_i$}, the prior $P(H_k)$ is the probability that a star belongs to hypothesis $H_k$ irrespective of its kinematic properties, and the fully marginalized likelihood

<footnote n="9">⁹ The Lupus star-forming region is slightly above this limit at ~155 pc (Lombardi et al. 2008), and the distance of the Chamaeleon star-forming region has recently been revised above 150 pc; (Voirin et al. 2017).</footnote>

<footnote n="10">¹⁰ Interactive Data Language.</footnote>

<footnote n="11">¹¹ IDL: https://github.com/jgagneastro/banyan_sigma_idl as archived in Zenodo (Gagné 2018b); Python: https://github.com/jgagneastro/banyan_sigma as archived in Zenodo (Gagné 2018a).</footnote>

<footnote n="12">¹² http://www.exoplanetes.umontreal.ca/banyan/banyansigma.php</footnote>
$P(\{O_i\})$ is the probability that a star displays observables $\{O_i\}$ irrespective of its membership. Once the priors and likelihoods are determined, the fully marginalized likelihood can be obtained with:

$$P(\{O_i\}) = \sum_k P(\{O_i\} | H_k) P(H_k). \quad (2)$$

It is often the case that the radial velocity ($\nu$) and/or distance ($\varpi$) of a star are not known, preventing a direct calculation of the likelihood $P(\{O_i\} | H_k)$. The case where both measurements are missing will be considered here. In this scenario, a likelihood probability can still be obtained by marginalizing over both missing parameters:

$$P(\{O_i\} | H_k) = \int_{-\infty}^{\infty} \int_{0}^{\infty} P_o(\{O_i\} | H_k) d\varpi \, d\nu, \quad (3)$$

where the symbol $P_o$ is used to distinguish probability densities from probabilities $P$ that are free of physical units.

Gagné et al. (2014) demonstrated that young associations can be well described with Gaussian distributions by working in the Galactic position XYZ and space velocity UVW frame of reference $\{Q_i\}$. The distributions of direct observables $\{O_i\}$ for the members of a young association would be accurately described only by complex functions in this coordinate frame. For this reason, BANYAN II approximated the likelihood by computing it directly in the $\{Q_i\}$ parameter space:

$$P(\{O_i\} | H_k) = \int_{-\infty}^{\infty} \int_{0}^{\infty} P_o(\{Q_j(\{O_i\}', \nu, \varpi)\} | H_k) d\varpi \, d\nu, \quad (4)$$

where $\{O_i\}'$ represents the set of observables excluding the radial velocity $\nu$ and distance $\varpi$. Equation (4) inherently ignores the Jacobian of the transformation $\{O_i\} \to \{Q_i\}$, discussed further in Section 3. In Gagné et al. (2014), both integrals were solved numerically on a uniform grid of 500 × 500 points over $\nu$ and $\varpi$, covering −35 to 35 km s⁻¹ and 0.1 to 200 pc, respectively. On each point ($\nu, \varpi$) of the grid, the observables $\{O_i\}'$ were transformed to the $\{Q_i\}$ frame of reference, and compared with a model of hypothesis $H_k$ to derive the probability density $P_o(\{Q_j(\{O_i\}', \nu, \varpi)\} | H_k)$. The sum of all probability densities on the grid were then taken as an approximation of the likelihood $P(\{O_i\} | H_k)$. These approximations were mostly limiting in that they prevented the inclusion of high-velocity ($|\nu| > 35$ km s⁻¹) or distant ($\varpi > 200$ pc) stars, but they also required 250,000 probability density calculations for each star.

As measurement errors on the sky position of a star are always small enough to have a negligible contribution to XYZUVW compared to those on proper motion, radial velocity or distance measurements, they were ignored completely in BANYAN II, and will still be ignored in BANYAN $\Sigma$. Properly including the measurement errors on proper motion would require the introduction of two additional integrals to obtain a modified likelihood that takes error bars into account:

$$P_e(\{O_i\} | H_k) = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} P(\{O_i\} | H_k) P_m(\mu_\alpha, \mu_\delta) d\mu_\alpha d\mu_\delta,$$

where $P_m(\mu_\alpha, \mu_\delta)$ is a probability density function describing the proper motion measurement, such as the product of two Gaussian distributions centered on the measured values with the appropriate characteristic widths. Because numerically

solving this likelihood would be impractical, BANYAN II approximated their effect by using a propagation error formula to obtain error bars on $U, V$, and $W$, and added them in quadrature to the characteristic widths of the Gaussian models describing each hypothesis $H_k$.

The kinematic models of BANYAN II are based on Gaussian ellipsoids in XYZ and UVW space, which are freely rotated along any axis. Models for seven young associations were included: TW Hya (TWA; de La Reza et al. 1989; Kastner et al. 1997), $\beta$ Pictoris ($\beta$PMG; Zuckerman et al. 2001a), Tucana-Horologium (THA; Torres et al. 2000; Zuckerman et al. 2001b), Carina (CAR; Torres et al. 2008), Columba (COL; Torres et al. 2008), Argus (ARG; Makarov & Urban 2000) and AB Doradus (ABDMG; Zuckerman et al. 2004). The model of field stars was obtained by fitting a spatial and a kinematic Gaussian ellipsoid to the Besançon Galactic model within 200 pc.

### 3. Banyan $\Sigma$: An Improved Algorithm

The framework of BANYAN $\Sigma$ improves on BANYAN II by (1) using the analytical solution to the marginalization integrals over radial velocity and distance, (2) using multivariate Gaussian models for the young associations and a mixture of multivariate Gaussians to model the Galactic field, (3) removing several approximations in the calculation of the Bayesian likelihood, (4) accounting for parallax motion, and (5) including a larger number of young associations. The algorithm of BANYAN $\Sigma$ is described in this section, and additional improvements with respect to the kinematic models are described in Sections 4–6. The models of BANYAN $\Sigma$ are built from a “training set” consisting of a set of bona fide or high-likelihood members of young associations compiled from the literature. The models are therefore not built statistically, and the algorithm of BANYAN $\Sigma$ is analogous to a Bayesian classification algorithm with a Gaussian mixture model (e.g., see Bishop & Nasrabadi 2007 and McLachlan & Peel 2000).

The multivariate Gaussian models of BANYAN $\Sigma$ are described in Section 3.1, and the change of coordinates from the direct observables $\{O_i\}$ to the Galactic frame of reference $\{Q_i\}$ is described in Section 3.2. The analytical solution of the Bayesian likelihood is presented in Section 3.3, and determination of the radial velocity $\nu$ and distance $\varpi$ that maximize the Bayesian likelihood (Equation (3)) are developed in Section 3.4. Section 3.5 presents a method to approximate the effect of measurement errors on proper motion. Section 3.6 details a new improvement to BANYAN $\Sigma$, allowing it to apply a parallax motion correction when the distance of a star is not known and its proper motion measurement is based on two epochs only. Sections 3.7–3.9 describe additional options to the BANYAN $\Sigma$ algorithm, to include measurements of radial velocity and/or distance, constraints from spectrophotometric observables, or to ignore the spatial distribution of young associations.

#### 3.1. Kinematic Models

The distribution of stars in young stellar associations and the Galactic field are modeled with multivariate Gaussian distributions in XYZUVW space. This is a generalization over the freely rotating individual XYZ and UVW Gaussian models used in both BANYAN II (Gagné et al. 2014) and LACEwING (Riedel et al. 2017b), as it models correlations between mixed spatial
and kinematic coordinates. The kinematic model corresponding to hypothesis $H_k$ can be written as:

$$ 
\begin{aligned}
\mathcal{P}_M(\bar{Q}, \bar{\tau}, \bar{\Sigma}) &= \frac{e^{-\frac{1}{2}\mathcal{M}^2}}{\sqrt{(2\pi)^6|\bar{\Sigma}|}}, \\
\text{where } \mathcal{M} &= \sqrt{(\bar{Q} - \bar{\tau})^T \bar{\Sigma}^{-1} (\bar{Q} - \bar{\tau})} \quad (5)
\end{aligned}
$$

is the Mahalanobis distance and the bar (e.g., $\bar{v}$) and double-bar (e.g., $\bar{M}$) symbols are used to indicate vectors and matrices, respectively, in XYZUVW space. $\bar{Q}$ is a 6D vector built from {$Q_i$}, which correspond to the XYZUVW coordinates of an object. The Mahalanobis distance is a generalization of the concept of measuring how many standard deviations a data point is from the center of a Gaussian distribution, and is applied to multivariate Gaussian distributions in the present work. A Mahalanobis distance has no units and accounts for correlations in the multivariate Gaussian probability density function (Mahalanobis 1936).

The multivariate Gaussian model includes a total of 27 free parameters: six are stored in the $\bar{\tau}$ vector and indicate the center of the association; six are stored in the diagonal of the covariance matrix $\bar{\Sigma}$ and indicate the 6D size of the association; and 15 more are stored in the independent off-diagonal elements of $\bar{\Sigma}$ and indicate the orientation of the ellipsoid in 6D space or, equivalently, the correlations between each combination of coordinates. $\bar{x}^T$ indicates a vector transposition and $\bar{s}^{-1}$ a matrix inverse.

The Bayesian likelihoods in the Galactic frame of reference {$Q_i$} can thus be written as:

$$ \mathcal{P}_q(\{Q_i\} | H_k) = \mathcal{P}_M(\bar{Q}, \bar{\tau}, \bar{\Sigma}), $$

where the dependencies of $\bar{\tau}$ and $\bar{\Sigma}$ on the association index $k$ are implicit.

A simple approach for obtaining the parameters of a kinematic model is to calculate the average position $\bar{\tau}$ of the members in XYZUVW space and their variances and covariances to build the $\bar{\Sigma}$ matrix. A method that is more robust to outliers, and accounts for individual measurement errors, is presented in Section 5.

## 3.2. Change of Coordinates

Solving the Bayesian likelihood in Equation (3) requires applying a change of coordinates from the observables frame of reference {$Q_i$} to the Galactic frame of reference {$Q_i$}. The equations for this transformation are detailed by Johnson & Soderblom (1987), where the components of $\bar{Q}_i$ can be written as a linear combination of the radial velocity $\nu$ and the distance $\varpi$:

$$ \bar{Q} = \bar{\Omega}\nu + \bar{\Gamma}\varpi; \quad (6) $$

the components of $\bar{\Omega}$ and $\bar{\Gamma}$ are:

$$ 
\begin{aligned}
\bar{\Omega} &= (0, 0, 0, M_0, M_1, M_2) = (0, \mathcal{M}), \\
\bar{\Gamma} &= (\lambda_0, \lambda_1, \lambda_2, N_0, N_1, N_2) = (\lambda, \mathcal{N}),
\end{aligned}
$$

and symbols in bold represent 3D vectors or matrices in XYZ or UVW space.

The vectors $\lambda$, $\mathcal{M}$ and $\mathcal{N}$ transform the sky position, proper motion, radial velocity, and distance to XYZUVW following:

$$ 
\begin{aligned}
\lambda &= (\cos b \cos l, \cos b \sin l, \sin b), \\
\mathcal{M} &= \mathcal{T}\mathcal{A}\mathcal{m}, \\
\mathcal{N} &= \mathcal{T}\mathcal{A}\mathcal{n}, \\
\mathcal{A} &= \begin{bmatrix} \cos \alpha \cos \delta & -\sin \alpha & -\cos \alpha \sin \delta \\ \sin \alpha \cos \delta & \cos \alpha & -\sin \alpha \sin \delta \\ \sin \delta & 0 & \cos \delta \end{bmatrix}, \\
\mathbf{m} &= (1, 0, 0), \\
\mathbf{n} &= \kappa(0, \mu_{\alpha} \cos \delta, \mu_{\delta}),
\end{aligned}
\quad (7) $$

where $l$ and $b$ are the Galactic longitude and latitude, $\alpha$ and $\delta$ are the R.A. and decl., $\mu_{\alpha} \cos \delta$ and $\mu_{\delta}$ are the proper motion, and $\kappa \approx 4.74 \cdot 10^{-3}$ corresponds to $10^{-3}$ au yr$^{-1}$ so that proper motions are expressed in mas yr$^{-1}$.

The $\mathcal{T}$ matrix is a combination of rotation matrices involving the equatorial position of the North Galactic Pole, and is detailed by Johnson & Soderblom (1987). We however use a definition of $\mathcal{T}$ where the first row has the opposite sign of that defined by Johnson & Soderblom (1987), so that $\mathcal{U}$ points toward the Galactic center and UVW forms a right-handed system:<fnref n="13" />

$$ \mathcal{T} = \begin{bmatrix} -0.054875560 & -0.87343709 & -0.48383502 \\ 0.49410943 & -0.44482963 & 0.74698224 \\ -0.86766615 & -0.19807637 & 0.45598378 \end{bmatrix}. $$

Equation (6) can be used to express $\mathcal{P}_q$ as a function of the observables {$O_i$}, but solving the marginalization integrals of Equation (3) requires applying a change of coordinates (from the {$Q_i$} frame to the {$O_i$} frame) to the probability density function itself $\mathcal{P}_q \to \mathcal{P}_o$. This step is detailed in Appendix A, and yields:

$$ \mathcal{P}_o(\{O_i\} | H_k) = \varpi^4 \mathcal{P}_q(\{Q_i\} | H_k). $$

Inserting the coordinate transformation (Equation (6)) into the kinematic model defined in Equation (5) yields:

$$ 
\begin{aligned}
\mathcal{M}_k^2(\nu, \varpi) = & \sum_{ij} [\sigma_{ij,k}^{-1} \Omega_i \Omega_j \nu^2 + \sigma_{ij,k}^{-1} \Gamma_i \Gamma_j \varpi^2 \\
& + \sigma_{ij,k}^{-1} (\Omega_i \Gamma_j + \Gamma_i \Omega_j) \nu \varpi \\
& - \sigma_{ij,k}^{-1} (\Omega_i \tau_{j,k} + \tau_{i,k} \Omega_i) \nu \\
& - \sigma_{ij,k}^{-1} (\Gamma_i \tau_{j,k} + \tau_{i,j} \Gamma_j) \varpi \\
& + \sigma_{ij,k}^{-1} \tau_{i,k} \tau_{j,k}].
\end{aligned}
\quad (8) $$

All terms in $\Omega_i$, $\Gamma_i$ and $\tau_{i,k}$ can be described as scalar products induced by the inverse covariance matrix $\sigma_{ij,k}^{-1}$, such as:

$$ \langle \bar{X}, \bar{Y} \rangle_k = \sum_{ij} \sigma_{ij,k}^{-1} X_i Y_j, \quad (9) $$

where the $k$ index will be omitted in the remainder of this work for simplicity.

Equation (8) can be further simplified from the fact that the covariance matrix is symmetric by definition. With these two

<footnote n="13">¹³ This display of the $\mathcal{T}$ matrix follows the column major convention (e.g., IDL or Fortran). Because Python follows the row major matrix convention, $\mathcal{T}$ must be transposed in that language. In IDL, the $\mathcal{T}$ matrix is encoded with $T = [[-0.0548\dots, 0.494\dots, \dots], [-0.873\dots, \dots], [\dots]]$. In Python, it is encoded with $T = np.array([[[-0.0548\dots, -0.873\dots, \dots], [0.494\dots, \dots], [\dots]]]).$</footnote>
simplifications we get:

$$ M^2(\nu, \bar{\omega}) = \langle \bar{\Omega}, \bar{\Omega} \rangle \nu^2 + \langle \bar{\Gamma}, \bar{\Gamma} \rangle \bar{\omega}^2 + 2 \langle \bar{\Omega}, \bar{\Gamma} \rangle \nu \bar{\omega} \\ - 2 \langle \bar{\Omega}, \bar{\tau} \rangle \nu - 2 \langle \bar{\Gamma}, \bar{\tau} \rangle \bar{\omega} + \langle \bar{\tau}, \bar{\tau} \rangle. \quad (10) $$

This completes the coordinate change of the Bayesian likelihood:

$$ P_o(\{O_i\}|H) d\nu d\bar{\omega} = \frac{\bar{\omega}^4 e^{-\frac{1}{2}M^2(\nu, \bar{\omega})}}{\sqrt{(2\pi)^6 |\bar{\Sigma}|}} d\nu d\bar{\omega}. \quad (11) $$

### 3.3. Solving the Marginalization Integrals

A complete analytical solution to Equation (3) is developed in Appendix B and yields:

$$ P(\{O_i\}|H) = \frac{D'_{-5}(\gamma/\sqrt{2\beta})e^{\gamma^2/4\beta - \zeta}}{|\bar{\Omega}|\sqrt{\pi^5\beta^5|\bar{\Sigma}|}}, \quad (12) $$

$$ \text{with } \beta = \frac{\langle \bar{\Gamma}, \bar{\Gamma} \rangle}{2} - \frac{1}{2} \frac{\langle \bar{\Omega}, \bar{\Gamma} \rangle^2}{\langle \bar{\Omega}, \bar{\Omega} \rangle}, \quad (13) $$

$$ \gamma = \frac{\langle \bar{\Omega}, \bar{\Gamma} \rangle \langle \bar{\Omega}, \bar{\tau} \rangle}{\langle \bar{\Omega}, \bar{\Omega} \rangle} - \langle \bar{\Gamma}, \bar{\tau} \rangle, \quad (14) $$

$$ \zeta = \frac{\langle \bar{\tau}, \bar{\tau} \rangle}{2} - \frac{1}{2} \frac{\langle \bar{\Omega}, \bar{\tau} \rangle^2}{\langle \bar{\Omega}, \bar{\Omega} \rangle}, \quad (15) $$

$$ D'_{-5}(x) = \sqrt{\frac{\pi}{2}}(x^4 + 6x^2 + 3) \operatorname{erfc}\left(\frac{x}{\sqrt{2}}\right) - (x^3 + 5x)e^{-x^2/2}, $$

where $D'_{-5}(x)$ is a parabolic cylinder function (Magnus & Oberhettinger 1948) modified for numerical stability.

A limitation of this development is that correlations between the measurement errors of the sky coordinates, proper motion and parallax cannot be accounted for. Properly accounting for such correlations would require performing much more CPU-intensive and less precise numerical integrals. However, we demonstrate in Section 8 that ignoring such correlations cause negligible biases in the Bayesian membership probabilities.

In summary, obtaining the Bayesian probability for a given star and hypothesis H requires calculating the components of the vectors $\bar{\Omega}$, $\bar{\Gamma}$, $\bar{\tau}$ (Equations (6) and (7)) and their various scalar products (Equations (13)–(15)) using Equation (9), then evaluating the non-marginalized Bayesian likelihood with Equation (12), and finally evaluating the Bayesian membership probability with Equation (1), which makes use of Equation (2). The priors for each hypothesis are also needed to evaluate Equation (12); these are defined in Section 7. When a large number of stars is analyzed, it is possible to improve the efficiency in solving Equation (12) by calculating the individual six components of the 6D vectors, as well as $\beta$, $\gamma$, $\zeta$ and $P(\{O_i\}|H)$ for the full array of stars at once.

### 3.4. Optimal Radial Velocity and Distance

The optimal values for the radial velocity and distance ($\nu_o$, $\bar{\omega}_o$) that maximize the non-marginalized Bayesian likelihood $P_o(\{O_i\}|H)$ can be determined for each hypothesis H. They correspond to predictions for the radial velocity and distance of a star, assuming that the star is a true member of H (e.g., Liu et al. 2016 obtain BANYAN II distance predictions with accuracies as low as ~20% by ignoring this fact).

The optimal radial velocity and distance are derived in Appendix C:

$$ \begin{aligned} \bar{\omega}_o &= \frac{-\gamma + \sqrt{\gamma^2 + 32\beta}}{4\beta}, \\ \nu_o &= \frac{4 + \langle \bar{\Gamma}, \bar{\tau} \rangle \bar{\omega}_o - \langle \bar{\Gamma}, \bar{\Gamma} \rangle \bar{\omega}_o^2}{\langle \bar{\Omega}, \bar{\Gamma} \rangle \bar{\omega}_o}, \\ \sigma_{\bar{\omega}} &= |\bar{\Gamma}|^{-1}, \\ \sigma_{\nu} &= |\bar{\Omega}|^{-1}, \end{aligned} $$

where $\sigma_{\bar{\omega}}$ and $\sigma_{\nu}$ represent statistical $1\sigma$ error bars on the optimal values. These two values are given for each hypothesis and each star as an output of BANYAN $\Sigma$, and can be taken as predictions on the radial velocity and distance measurements that would maximize the probability of a given hypothesis.

Because the optimal radial velocity and distance are intrinsically linked to the assumption that a star is a member of hypothesis H, adopting the values ($\nu_o$, $\bar{\omega}_o$) for a low-probability hypothesis H inherently carries a large risk of the true radial velocity and distance of a star of being several $\sigma$ away from the prediction. The values of $\sigma_{\bar{\omega}}$ and $\sigma_{\nu}$ will therefore be unreliable in such a situation.

### 3.5. Approximating the Effect of Proper Motion Measurement Errors

As discussed in Section 2, including the effect of proper motion errors requires solving two additional marginalization integrals on the proper motion components, in addition to those on radial velocity and distance. Since obtaining an analytical solution to these four marginalization integrals is impractical, this section describes an analytical approximation for the effect of the proper motion measurement error.

The optimal radial velocity $\nu_o$ and distance $\bar{\omega}_o$ that maximize the non-marginalized Bayesian likelihood $P_o(\{O_i\}|H)$ can be used to propagate the proper motion measurement errors to the galactic position XYZ and space velocity UVW in the vicinity of the maximum of the membership probability distribution function. The sky position, proper motion, $\nu_o$, $\bar{\omega}_o$ and proper motion errors are propagated in XYZUVW space to obtain a Galactic error vector $\bar{\delta}_q = (\sigma_X, \sigma_Y, \sigma_Z, \sigma_U, \sigma_V, \sigma_W)$.

It is then possible to add these errors in quadrature to the diagonal elements of $\bar{\Sigma}$ without affecting the orientation of the multivariate Gaussian with:

$$ \begin{aligned} \bar{\Sigma}' &= \bar{\bar{G}} \bar{\Sigma} \bar{\bar{G}}, \\ \bar{\bar{G}} &= \operatorname{diag}_+ \sqrt{\frac{\operatorname{diag}_-(\bar{\Sigma}) + \bar{\sigma}_q^2}{\operatorname{diag}_-(\bar{\Sigma})}}, \end{aligned} $$

where $\operatorname{diag}_-(\bar{\bar{M}})$ extracts the diagonal elements of matrix $\bar{\bar{M}}$ and $\operatorname{diag}_+(\bar{\bar{v}})$ builds a diagonal matrix with vector $\bar{\bar{v}}$.

The need for a time-consuming matrix inversion of $\bar{\Sigma}'$ for each star to evaluate Equation (10) can be avoided with:

$$ \begin{aligned} \bar{\Sigma}'^{-1} &= \bar{\bar{G}}^{-1} \bar{\Sigma}'^{-1} \bar{\bar{G}}^{-1}, \\ \bar{\bar{G}}^{-1} &= \operatorname{diag}_+ \sqrt{\frac{\operatorname{diag}_-(\bar{\Sigma})}{\operatorname{diag}_-(\bar{\Sigma}) + \bar{\sigma}_q^2}}. \end{aligned} $$
Including an approximated effect of the proper motion measurement errors will thus require a calculation of (1) the $\bar{\Omega}$, $\bar{\Gamma}$ and $\bar{\tau}$ vectors and their scalar products, (2) $\nu_o$ and $\bar{\omega}_o$, (3) the $\bar{\sigma}_q$ vector and the corresponding inflated matrix $\bar{\Sigma}'$, and (4) all quantities from the beginning of Section 3.2 obtained with the updated scalar product based on $\bar{\Sigma}'^{-1}$. The approximation described here does not make the assumption that a star is a member of any young association or the field. Instead, the proper motion errors are propagated to XYZUVW independently for each Bayesian hypothesis, ensuring that it is valid near the peak of the probability distribution of each hypothesis. In other words, the steps described above are carried out independently when calculating the membership probability of each hypothesis.

### 3.6. Parallax Motion

When the proper motion of a nearby star is measured based on two epochs only, the measurement may be contaminated in part by parallax motion. This is true because the measurement of a star's displacement between two epochs include the compounding effects of its true proper motion (i.e., its space velocity projected on the celestial sphere) with that of its displacement along its parallactic ellipse, the latter of which is caused by a change in the observer's point of view as the Earth progresses on its orbit around the Sun. These two effects can however be decoupled in the BANYAN $\Sigma$ formalism, as long as the parallax factors ($\psi_\alpha, \psi_\delta$), described below, are measured for the star. The parallax factors physically represent the motion of the star purely due to the Earth's motion between the two epochs if the star was placed at exactly 1 pc from the Sun.

The parallax motion ($\Delta_{\alpha\pi}, \Delta_{\delta\pi}$) of a star is given by Smart & Green (1977):

$$ \Delta_{\alpha\pi} \cos \delta = \frac{\cos \alpha \cos e \sin \ell_s - \sin \alpha \cos \ell_s}{\bar{\omega}}, \quad (16) $$

$$ \Delta_{\delta\pi} = \frac{\cos \delta \sin e \sin \ell_s - \cos \alpha \sin \delta \cos \ell_s}{\bar{\omega}} - \frac{\sin \alpha \sin \delta \cos e \sin \ell_s}{\bar{\omega}}, \quad (17) $$

where $e$ is the obliquity of the ecliptic of the Earth's orbit and $\ell_s$ is the ecliptic longitude of the Sun at a given epoch.<fnref n="14" /> These equations can be simplified by grouping all epoch-dependent terms into ($\phi_\alpha, \phi_\delta$):

$$ \Delta_{\alpha\pi}(\alpha, t) \cos \delta = \frac{\phi_{\alpha}(\alpha, t)}{\bar{\omega}}, \quad (18) $$

$$ \Delta_{\delta\pi}(\alpha, \delta, t) = \frac{\phi_{\delta}(\alpha, \delta, t)}{\bar{\omega}}, \quad (19) $$

where $t$ is the epoch.

The apparent motion ($\mu'_\alpha, \mu'_\delta$) of a star between epochs $t_1$ and $t_2$ will thus be given by:

$$ \begin{aligned} \mu'_{\delta} &= \frac{(\delta(t_2) + \Delta_{\delta\pi}(t_2)) - (\delta(t_1) + \Delta_{\delta\pi}(t_1))}{t_2 - t_1} \\ &= \mu_{\delta} + \frac{1}{\bar{\omega}} \frac{\phi_{\delta}(t_2) - \phi_{\delta}(t_1)}{t_2 - t_1} \\ &= \mu_{\delta} + \frac{\psi_{\delta}}{\bar{\omega}}, \quad \text{and similarly:} \end{aligned} $$

$$ \mu'_{\alpha} \cos \delta = \mu_{\alpha} \cos \delta + \frac{\psi_{\alpha}}{\bar{\omega}}. $$

Since Equation (6) is linear in proper motion components, it can be expressed as a function of apparent motion with the form:

$$ \begin{aligned} \bar{Q} &= \bar{\Omega}\nu + \bar{\Gamma}'\bar{\omega} - \bar{\Phi}, \\ \bar{\Phi} &= (0, \mathcal{T}\mathcal{A}\psi), \\ \psi &= \kappa(0, \psi_\alpha, \psi_\delta), \\ \psi_\alpha &= \frac{\bar{\omega} \cos \delta (\Delta_{\alpha\pi}(t_2) - \Delta_{\alpha\pi}(t_1))}{t_2 - t_1}, \\ \psi_\delta &= \frac{\bar{\omega} (\Delta_{\delta\pi}(t_2) - \Delta_{\delta\pi}(t_1))}{t_2 - t_1}, \end{aligned} \quad (20) $$

where $\bar{\Gamma}'$ is a function of apparent motion ($\mu'_\alpha \cos \delta, \mu'_\delta$), the quantity that is directly measured, instead of true proper motion ($\mu_\alpha \cos \delta, \mu_\delta$).

Because $\bar{Q}$ only appears in the Bayesian likelihood as relative to the center of the moving group model $\tau$, the effect of parallax motion can be fully accounted for by shifting the UVW center of the young association kinematic model by $+\mathcal{T}\mathcal{A}\psi$, which is equivalent to shifting $\bar{\tau}$ by $+\bar{\Phi}$:

$$ \bar{\tau}' = \bar{\tau} + \bar{\Phi}. $$

As a consequence, the parallax motion can be accounted for by using the measured apparent motion as if it were a true proper motion, and replacing and $\bar{\tau} \to \bar{\tau}'$ in the BANYAN $\Sigma$ formalism. This requires measurements of the parallax factors ($\psi_\alpha, \psi_\delta$) for each star in addition to measurements of their apparent motion. In practice, this correction is applied by the BANYAN $\Sigma$ software only when the keyword `use_psi` is explicitly used. This indicates that: (1) the proper motion that is input to BANYAN $\Sigma$ was measured from two epochs only, (2) the effect of parallax motion was not corrected in the proper motion measurement and therefore it really is a measurement of apparent motion, and (3) the parallax factors are provided to BANYAN $\Sigma$ using the same two epochs as those between which the proper motion was measured. If any of the above statements are not true, the parallax motion correction described above should not be used.

### 3.7. Additional Kinematic Observables

Only sky position and proper motion are required for BANYAN $\Sigma$ to compute membership probabilities. However, radial velocities and/or distances can also be included as input measurements to obtain more accurate membership probabilities. This is similar to the functioning of BANYAN I (Malo et al. 2013) and BANYAN II (Gagné et al. 2014). In cases where a radial velocity measurement is available, Equation (3)

<footnote n="14">¹⁴ $e$ and $\ell_s$ can be computed from the sky position $\alpha, \delta$ and the julian date using the IDL astrolib routine `sunpos.pro`.</footnote>
can be rewritten as:

$$P(\{O_i\}|H) = \int_{-\infty}^{\infty} P_m(\nu) \int_{0}^{\infty} \frac{\varpi^4 e^{-\frac{1}{2}\mathcal{M}^2}}{\sqrt{(2\pi)^6 |\Sigma'|}} d\varpi d\nu,$$

where $P_m(\nu)$ is the probability density function that represents the radial velocity measurement.

Assuming that $P_m(\nu)$ is a Gaussian distribution centered on $\nu_m$ with a characteristic width of $\sigma_{\nu,m}$, the equation above can be solved with a similar method to that described in Appendix B, where the following scalar products are modified with:

$$\begin{aligned} \langle \bar{\Omega}, \bar{\Omega} \rangle &\rightarrow \langle \bar{\Omega}, \bar{\Omega} \rangle + (\sigma_{\nu,m})^{-2}, \\ \langle \bar{\Omega}, \bar{\tau} \rangle &\rightarrow \langle \bar{\Omega}, \bar{\tau} \rangle + \nu_m (\sigma_{\nu,m})^{-2}, \\ \langle \bar{\tau}, \bar{\tau} \rangle &\rightarrow \langle \bar{\tau}, \bar{\tau} \rangle + \nu_m^2 (\sigma_{\nu,m})^{-2}. \end{aligned}$$

The case where a distance measurement is available can be solved in a similar way, by replacing:

$$\begin{aligned} \langle \bar{\Gamma}, \bar{\Gamma} \rangle &\rightarrow \langle \bar{\Gamma}, \bar{\Gamma} \rangle + (\sigma_{\varpi,m})^{-2}, \\ \langle \bar{\Gamma}, \bar{\tau} \rangle &\rightarrow \langle \bar{\Gamma}, \bar{\tau} \rangle + \varpi_m (\sigma_{\varpi,m})^{-2}, \\ \langle \bar{\tau}, \bar{\tau} \rangle &\rightarrow \langle \bar{\tau}, \bar{\tau} \rangle + \varpi_m^2 (\sigma_{\varpi,m})^{-2}. \end{aligned}$$

The case where both radial velocity and distance measurements are available can be solved by combining all of the variable changes described above (the two changes on $\langle \bar{\tau}, \bar{\tau} \rangle$ must be cumulated).

### 3.8. Photometric Observables

It is possible to constrain the distance of a star from its position in a color-magnitude or spectral type-magnitude diagram by comparing its absolute magnitude to a sequence of field stars, or to members of a young association, at a fixed color or spectral type. The position of a sequence in most of these diagrams is dependent on the age of its population, which translates to a different distance constraint for each Bayesian hypothesis.

Such photometric constraints can be included in the BANYAN $\Sigma$ framework, in a similar way to the method described in Section 3.7 for distance measurements, except that different values of the most likely distance $\varpi_m$ and its uncertainty $\sigma_{\varpi,m}$ must be used, one for each hypothesis, because they derive from different color-magnitude sequences. This is a consequence of the fact that the different Bayesian hypotheses correspond to populations of stars at different ages.

In the absence of a trigonometric distance measurement, users can create custom color-magnitude diagrams and determine values of $(\varpi_m, \sigma_{\varpi,m})$ for the field and each young association, and provide them to BANYAN $\Sigma$ for a full inclusion of these constraint in the Bayesian probabilities. Multiple color-magnitude diagrams can also be combined into single measurements of $(\varpi_m, \sigma_{\varpi,m})$ for a given star and young association, but failing to account for covariances between different photometric bands in a given stellar population would result in artificially small values of $\sigma_{\varpi,m}$. Such unrealistically precise constraints on the distance of a star would hinder the ability of BANYAN $\Sigma$ to correctly identify the candidate members of a young association. Both the IDL and Python implementations of BANYAN $\Sigma$ can accept these photometric distance constraints through the `constraint_dist_per_hyp` and `constraint_edist_per_hyp` keywords, which are detailed in the documentation of the code.

In summary, distinct color-magnitude diagrams for each hypothesis can be included by providing BANYAN $\Sigma$ with distinct photometric constraints on the distance of a star.

In the cases where both parallax and photometric measurements are available, they can be included in a more straightforward way to BANYAN $\Sigma$ through the Bayesian priors: the vertical distances in a color-magnitude diagram between the measured absolute magnitude and the sequence of field or young objects at a fixed color can be transformed to a probability for each association using Bayes' theorem, and the natural logarithm of these photometric probabilities can be included in BANYAN $\Sigma$ with the `ln_priors` keyword available in both the IDL and Python implementations of the code.

Other age-dating observables, such as X-ray, UV, H$\alpha$, rotation and lithium abundance measurements can similarly be translated to a membership probability at the age of each young association (e.g., by comparing measurements with the X-ray luminosity distributions of Malo et al. 2014), and can also be included in the BANYAN $\Sigma$ prior probabilities. This framework allows users to add observables in the BANYAN $\Sigma$ membership determination without needing to change its algorithm, and remains accurate as long as no kinematic measurements are used to assign prior probabilities. The inclusion of such age indicators must rely on a user-specified method to translate each measurements into a probability that a given star is a member of each hypothesis, given the age of each young association. A compounded probability that each star is a member of each young association must then be calculated, based on only these youth indicators and no kinematics (e.g., by multiplying together the probabilities obtained from independent age indicators). The natural logarithm of these probabilities must then be input to BANYAN $\Sigma$ with the keyword `ln_priors`. These data are passed to BANYAN $\Sigma$ using a Python dictionary or an IDL structure depending on which version of the code is used, and we refer the reader to the respective documentation, which are provided as additional material to this manuscript, for more detail.

No color-magnitude sequences are provided here with the first version of BANYAN $\Sigma$, but they will be provided in future work as they are developed to target specific types of members.

### 3.9. Ignoring the Galactic Position XYZ

It is possible that the full spatial extent of some young associations have not yet been completely explored, especially at larger distances not covered by the Hipparcos survey. This possibility has been hypothesized by Bowler et al. (2017) among others, and the fact that the BANYAN tools rely on XYZ as well as UVW prevents an exploration of that potentially missing population of members.

The BANYAN $\Sigma$ framework can be adapted to rely on UVW only, by artificially setting the first three diagonal elements of the covariance matrices $\Sigma$ to a large value (e.g., $10^9$ pc$^2$), and setting all other elements that contain at least one spatial component to zero. This approach is equivalent to using very large spatial widths for all models of young associations and the field in the BANYAN $\Sigma$ formalism, and the general solution presented in Equation (12) remains unchanged. An option is provided in BANYAN $\Sigma$ to ignore the spatial XYZ coordinates, and can be used to locate young association
members that are spatially distant to the locus of known
members. However, the rate of contamination from field stars is
~100 times larger when using this option (see Section 8), and
we therefore recommend extreme caution when using it.

4. Bona Fide Members of Young Associations within 150 pc

In this section, a list of bona fide members of young
associations within 150 pc is compiled. This list will constitute
the training set for the Gaussian models used in BANYAN Σ (see
Sections 3.1 and 5 for more detail on the models). Each
association considered in this work is listed in Table 1 with its age
estimate and the total number of bona fide members that were
compiled. In the literature, stars are typically considered bona fide
members when they benefit from signs of youth and full
kinematic measurements that allow them to be placed in XYZUVW
space (Malo et al. 2013; Gagné et al. 2014). The indicators of
youth used in the literature depend on the spectral type of the
stars, and include isochronal age determinations through color-
magnitude positions, lithium measurements, UV or X-ray
luminosity, rotational velocity, Hα emission, and rotational
velocity; see Soderblom et al. (2014) for a review of these age-
dating methods. Here we require the same measurements for bona
fide members of the nearest or most well-studied young
associations. Nine of the 27 associations that are further away
than ~90 pc do not have enough members with full 6D
kinematics to require them in the construction of their model. In
these cases, an average radial velocity or parallax (or both) for the
association are used instead of individual measurements.

In Section 4.1, we provide a summarized description of the 18 young associations for which models are built using the individual 6D kinematics of all members. The nine associations with incomplete kinematics are described in Section 4.2 and their adopted average distances and radial velocities are listed in Table 2. A description the Argus association, included in BANYAN II but excluded here, is provided in Section 4.3. A few individual objects that require further attention are discussed in Section 4.4. The addition of Gaia-DR1 data to several previously recognized high-probability candidate members of the young associations studied here makes them new bona fide members: a description of these new members is provided in Section 4.5. The method that we used to calculate 6D kinematics and their error bars from kinematic observables is described in Section 4.6.

All members compiled in this section were cross-matched
with the 2MASS, AllWISE (Wright et al. 2010; Kirkpatrick
et al. 2014) and the Gaia-DR1 catalogs. When available, sky
positions, proper motions and parallaxes from the Gaia-DR1
catalog were preferred to literature measurements. Targets
with no radial velocity measurements reported in the literature
membership lists were cross-matched with various catalogs
that provide radial velocities (Upgren & Harlow 1996;
Hawley et al. 1997; Bobylev et al. 2006; Torres et al. 2006;
Kharchenko et al. 2007; White et al. 2007; Fernández
et al. 2008; Shkolnik et al. 2011; Chubak et al. 2012;
Kordopatis et al. 2013; Malo et al. 2014). When authors did
not report radial velocity measurement errors, we adopted
those calculated by Riedel et al. (2017b; see their Table 6).

4.1. Associations With Full Kinematics

In this section, we provide a short description of the 18
young associations included in BANYAN Σ for which the

models will be built only from their *classical* bona fide
members, in the sense that only members with signs of youth
and full 6D kinematics are included.

The list of bona fide members presented in Gagné et al. (2014), which was largely based on that of Malo et al. (2013), was used as a starting point in this work. The Gagné et al. (2014) list includes members of TWA, βPMG, THA, CAR, COL, ARG and ABDMG. Gagné et al. (2017b) compiled an updated list of TWA members with new available data, and rejected contaminants from more distant associations, and was therefore adopted here. We refer the reader to Zuckerman & Song (2004) and Torres et al. (2008) for a detailed description of these associations.

Carina-Near (CARN) has been identified as a co-moving
group of ~200 Myr old stars by Zuckerman et al. (2006),
which includes a core of eight members and ten members that
are part of a spatially larger stream. Few studies have focused
on this moving group since its discovery, likely because it is
among the older ones. Gagné et al. (2017a) used a preliminary
version of BANYAN Σ to show that the variable T2.5 brown
dwarf SIMP J013656.5+093347 is likely a ~13 M<sub>jup</sub> member
of the CARN stream, and Riedel et al. (2017b) included CARN
in a young association classification tool (LACEwING) for the
first time.

The Ursa Major cluster (UMA; e.g., Eggen 1992) is a well-studied population of co-moving stars, consisting of a core of coeval young stars, and a stream of stars with heterogeneous compositions and ages. Soderblom & Mayor (1993) estimated an age of ~300 Myr for the core population, and Jones et al. (2015a, 2017) estimated an age of 414 ± 23 Myr based on interferometric measurements of its A-type members. The core membership lists of King et al. (2003) was adopted for this work, and the stream was not included in BANYAN Σ because of its heterogeneous nature.

The Hyades cluster (HYA) is a nearby (40–50 pc) and relatively young (600–800 Myr; Perryman et al. 1998) cluster that has been extensively studied in the literature (e.g., Perryman et al. 1998; Zuckerman & Song 2004). The membership list of Perryman et al. (1998) was adopted here.

The Upper Scorpius, (USCO), Upper Centaurus-Lupus (UCL) and Lower Centaurus-Crux (LCC) groups are part of the Sco-Cen star-forming region (Blaauw 1946; de Zeeuw et al. 1999), which consists of 5–30 Myr stars located at distances of ~110–150 pc. The membership lists of Rizzuto et al. (2011), Pecaut & Mamajek (2016) and Donaldson et al. (2017) were adopted here. Rizzuto et al. (2011) only list membership probabilities for the Sco-Cen region, and do not classify their members in the three subgroups. Their list was therefore cross-matched with that of de Zeeuw et al. (1999) to assign the correct sub-group, but all members of Sco-Cen that were newly discovered by Rizzuto et al. (2011) were not included at this stage. Once completed, the BANYAN Σ tool can be used to assign these new members to the correct subgroup; this is done in Section 9. Several radial velocities for USCO, UCL and LCC cataloged by Kharchenko et al. (2007) —which seem to be mistakenly listed as originating from Gontcharov (2006)—are astrometric<fnref n="15" /> radial velocities assuming moving group membership and an average UVW velocity. These measurements were rejected from our compilation.

<footnote n="15">¹⁵ Listed in the Kharchenko et al. (2007) catalog as “Index of Radial Velocity Catalogs” = 2 in table III/254/crvad2 of VizieR.</footnote>
<table><caption>**Table 1**

General Characteristics and Bayesian Priors of Young Associations</caption><thead><tr><th rowspan="2">Asso.</th><th rowspan="2">N<sub>k</sub><sup>a</sup></th><th colspan="5">ln α<sub>k</sub><sup>b</sup></th><th rowspan="2">⟨ω⟩<sup>c</sup><br>(pc)</th><th rowspan="2">⟨ν⟩<sup>d</sup><br>(km s<sup>-1</sup>)</th><th rowspan="2">S<sub>spa</sub><sup>e</sup><br>(pc)</th><th rowspan="2">S<sub>kin</sub><sup>f</sup><br>(km s<sup>-1</sup>)</th><th rowspan="2">Age<br>(Myr)</th><th rowspan="2">Age<br>Ref.</th></tr><tr><th>μ</th><th>μ, ν</th><th>μ, ω</th><th>μ, ν, ω</th></tr></thead><tbody><tr><td>118TAU</td><td>10</td><td>-17.22</td><td>-18.60</td><td>-21.37</td><td>-22.66</td><td>100 ± 10</td><td>14 ± 2</td><td>3.4</td><td>2.1</td><td>~10</td><td>1</td></tr><tr><td>ABDMG</td><td>48</td><td>-14.11</td><td>-15.39</td><td>-16.56</td><td>-17.60</td><td>30<sup>+20</sup><sub>-10</sub></td><td>10<sup>+10</sup><sub>-20</sub></td><td>19.0</td><td>1.4</td><td>149<sup>+51</sup><sub>-19</sub></td><td>2</td></tr><tr><td>βPMG</td><td>42</td><td>-13.57</td><td>-14.77</td><td>-17.39</td><td>-18.24</td><td>30<sup>+20</sup><sub>-10</sub></td><td>10 ± 10</td><td>14.8</td><td>1.4</td><td>24 ± 3</td><td>2</td></tr><tr><td>CAR</td><td>7</td><td>-13.41</td><td>-14.82</td><td>-18.45</td><td>-19.15</td><td>60 ± 20</td><td>20 ± 2</td><td>11.8</td><td>0.8</td><td>45<sup>+11</sup><sub>-7</sub></td><td>2</td></tr><tr><td>CARN</td><td>13</td><td>-15.51</td><td>-16.85</td><td>-17.64</td><td>-18.55</td><td>30 ± 20</td><td>15<sup>+7</sup><sub>-10</sub></td><td>14.0</td><td>2.1</td><td>~200</td><td>3</td></tr><tr><td>CBER</td><td>40</td><td>-13.70</td><td>-15.09</td><td>-22.32</td><td>-23.43</td><td>85<sup>+4</sup><sub>-5</sub></td><td>-0.1 ± 0.8</td><td>3.6</td><td>0.5</td><td>562<sup>+98</sup><sub>-84</sub></td><td>4</td></tr><tr><td>COL</td><td>23</td><td>-13.08</td><td>-14.10</td><td>-17.74</td><td>-18.34</td><td>50 ± 20</td><td>21<sup>+3</sup><sub>-8</sub></td><td>15.8</td><td>0.9</td><td>42<sup>+6</sup><sub>-4</sub></td><td>2</td></tr><tr><td>CRA</td><td>12</td><td>-17.55</td><td>-19.07</td><td>-21.89</td><td>-22.89</td><td>139 ± 4</td><td>-1 ± 1</td><td>1.5</td><td>1.7</td><td>4–5</td><td>5</td></tr><tr><td>EPSC</td><td>25</td><td>-17.47</td><td>-18.59</td><td>-22.38</td><td>-22.79</td><td>102 ± 4</td><td>14 ± 3</td><td>2.8</td><td>1.8</td><td>3.7<sup>+4.6</sup><sub>-1.4</sub></td><td>6</td></tr><tr><td>ETAC</td><td>16</td><td>-20.19</td><td>-21.36</td><td>-25.75</td><td>-26.22</td><td>95 ± 1</td><td>20 ± 3</td><td>0.6</td><td>2.0</td><td>11 ± 3</td><td>2</td></tr><tr><td>HYA</td><td>177</td><td>-20.02</td><td>-21.54</td><td>-22.14</td><td>-23.57</td><td>42 ± 7</td><td>39<sup>+3</sup><sub>-4</sub></td><td>4.5</td><td>1.2</td><td>750 ± 100</td><td>7</td></tr><tr><td>IC2391</td><td>16</td><td>-18.05</td><td>-18.90</td><td>-21.55</td><td>-21.55</td><td>149 ± 6</td><td>15 ± 3</td><td>2.2</td><td>1.4</td><td>50 ± 5</td><td>8</td></tr><tr><td>IC2602</td><td>17</td><td>-15.33</td><td>-16.60</td><td>-22.50</td><td>-22.60</td><td>146 ± 5</td><td>17 ± 3</td><td>1.8</td><td>1.1</td><td>46<sup>+6</sup><sub>-5</sub></td><td>9</td></tr><tr><td>LCC</td><td>82</td><td>-13.13</td><td>-14.27</td><td>-17.76</td><td>-17.77</td><td>110 ± 10</td><td>14 ± 5</td><td>11.6</td><td>2.2</td><td>15 ± 3</td><td>10</td></tr><tr><td>OCT</td><td>14</td><td>-11.74</td><td>-11.56</td><td>-13.85</td><td>-11.29</td><td>130<sup>+30</sup><sub>-20</sub></td><td>8<sup>+8</sup><sub>-9</sub></td><td>22.4</td><td>1.3</td><td>35 ± 5</td><td>11</td></tr><tr><td>PL8</td><td>11</td><td>-13.28</td><td>-14.54</td><td>-19.37</td><td>-19.75</td><td>130 ± 10</td><td>22 ± 2</td><td>5.0</td><td>1.1</td><td>~60</td><td>12</td></tr><tr><td>PLE</td><td>190</td><td>-18.72</td><td>-20.06</td><td>-20.61</td><td>-21.46</td><td>134 ± 9</td><td>6 ± 2</td><td>4.1</td><td>1.4</td><td>112 ± 5</td><td>13</td></tr><tr><td>ROPH</td><td>186</td><td>-17.49</td><td>-19.04</td><td>-24.10</td><td>-25.55</td><td>131 ± 1</td><td>-6.3 ± 0.2</td><td>0.7</td><td>1.6</td><td>&lt; 2</td><td>14</td></tr><tr><td>TAU</td><td>122</td><td>-10.39</td><td>-11.37</td><td>-17.04</td><td>-17.99</td><td>120 ± 10</td><td>16 ± 3</td><td>10.7</td><td>3.6</td><td>1–2</td><td>15</td></tr><tr><td>THA</td><td>39</td><td>-16.48</td><td>-17.78</td><td>-19.58</td><td>-20.25</td><td>46<sup>+8</sup><sub>-6</sub></td><td>9<sup>+5</sup><sub>-6</sub></td><td>9.1</td><td>0.8</td><td>45 ± 4</td><td>2</td></tr><tr><td>THOR</td><td>35</td><td>-14.05</td><td>-15.57</td><td>-20.22</td><td>-21.12</td><td>96 ± 2</td><td>19 ± 3</td><td>3.9</td><td>2.1</td><td>22<sup>+4</sup><sub>-3</sub></td><td>2</td></tr><tr><td>TWA</td><td>23</td><td>-16.99</td><td>-18.22</td><td>-20.42</td><td>-20.93</td><td>60 ± 10</td><td>10 ± 3</td><td>6.6</td><td>1.5</td><td>10 ± 3</td><td>2</td></tr><tr><td>UCL</td><td>103</td><td>-11.70</td><td>-13.10</td><td>-15.91</td><td>-16.19</td><td>130 ± 20</td><td>5 ± 5</td><td>17.4</td><td>2.5</td><td>16 ± 2</td><td>10</td></tr><tr><td>UCRA</td><td>10</td><td>-15.85</td><td>-16.87</td><td>-20.24</td><td>-20.48</td><td>147 ± 7</td><td>-1 ± 3</td><td>4.5</td><td>1.8</td><td>~10</td><td>16</td></tr><tr><td>UMA</td><td>9</td><td>-23.14</td><td>-24.01</td><td>-26.44</td><td>-27.13</td><td>25.4<sup>+0.8</sup><sub>-0.7</sub></td><td>-12 ± 3</td><td>1.2</td><td>1.3</td><td>414 ± 23</td><td>17</td></tr><tr><td>USCO</td><td>84</td><td>-12.77</td><td>-13.71</td><td>-17.62</td><td>-17.96</td><td>130 ± 20</td><td>-5 ± 4</td><td>9.9</td><td>2.8</td><td>10 ± 3</td><td>10</td></tr><tr><td>XFOR</td><td>11</td><td>-19.37</td><td>-20.80</td><td>-23.43</td><td>-23.72</td><td>100 ± 6</td><td>19 ± 2</td><td>2.6</td><td>1.3</td><td>~500</td><td>18</td></tr></tbody></table>

Notes. See Sections 4 and 5 for more detail. The full names of young associations are: 118 Tau (118TAU), AB Doradus (ABDMG), β Pictoris (βPMG), Carina (CAR), Carina-Near (CARN), Coma Berenices (CBER), Columba (COL), Corona Australis (CRA), ϵ Chamaeleontis (EPSC), η Chamaeleontis (ETAC), the Hyades cluster (HYA), Lower Centaurus Crux (LCC), Octans (OCT), Platais 8 (PL8), the Pleiades cluster (PLE), ρ Ophiuchi (ROPH), the Tucana-Horologium association (THA), 32 Orionis (THOR), TW Hya (TWA), Upper Centaurus Lupus (UCL), Upper CrA (UCRA), the core of the Ursa Major cluster (UMA), Upper Scorpius (USCO), Taurus (TAU) and χ<sup>1</sup> For (XFOR).

<sup>a</sup> Number of bona fide members included in the kinematic model.

<sup>b</sup> Bayesian prior ensuring a recovery rate of ≈50%–90% for a threshold *P* = 90%, depending on input observables.

<sup>c</sup> Peak of distance distribution and ±1σ range.

<sup>d</sup> Peak of radial velocity distribution and ±1σ range.

<sup>e</sup> Characteristic spatial scale in XYZ space.

<sup>f</sup> Characteristic kinematic scale in UVW space.

**References.** (1) Mamajek (2016), (2) Bell et al. (2015), (3) Zuckerman et al. (2006), (4) Silaj & Landstreet (2014), (5) Gennaro et al. (2012), (6) Murphy et al. (2013), (7) Brandt & Huang (2015), (8) Barrado y Navascués et al. (2004), (9) Dobbie et al. (2010), (10) Pecaut & Mamajek (2016), (11) Murphy & Lawson (2015), (12) Platais et al. (1998), (13) Dahm (2015), (14) Wilking et al. (2008), (15) Kenyon & Hartmann (1995), (16) this paper, (17) Jones et al. (2015b), (18) Pöhnl & Paunzen (2010).

The Octans association (OCT; Torres et al. 2008) is a group of young stars at ≈120 pc from the Sun, which has not been characterized as well as other young associations mainly due to its sky position located far in the Southern hemisphere (decl. between −87 and −20°). Murphy & Lawson (2015) performed a survey of its low-mass members and determined a lithium age of 20–40 Myr for this group. The members of OCT were compiled from Murphy & Lawson (2015).

The Pleiades cluster (PLE; Cummings 1921; Stauffer et al. 1989) is one of the best-studied clusters in the solar neighborhood. It is located at a distance of ~130 pc and recent estimates of its age based on its lithium depletion boundary are in the range ~110–120 Myr (Dahm 2015). The membership lists of Stauffer et al. (2007) and Galli et al. (2017) were adopted

here. Sarro et al. (2014) presented a Bayesian method to identify members of the pleiades based on multivariate Gaussians mixture models. This method differs from BANYAN Σ in that it does not consider other young associations, includes various photometric colors, works directly in proper motion space, and does not consider sky position because they study stars in the direction of the cluster only. The larger number of free parameters introduced by a mixture of multivariate Gaussians makes it possible to model the proper motion and color–magnitude distribution of the Pleiades members, which are not well represented by a single Gaussian distribution. The large number of known Pleiades members allows such a highly parametrized model, but it would likely be challenging to apply this method to sparser or nearby young associations.
<table><caption>**Table 2**

Adopted Average Distances and Radial Velocities for Young Associations with Partial Kinematics</caption><thead><tr><th>Association</th><th>ν<sub>avg</sub><br>(km s<sup>-1</sup>)</th><th>σ<sub>ν</sub><br>(km s<sup>-1</sup>)</th><th>N<sub>ν</sub></th><th>ω<sub>avg</sub><br>(pc)</th><th>σ<sub>ω</sub><br>(pc)</th><th>N<sub>ω</sub></th></tr></thead><tbody><tr><td>EPSC</td><td>...</td><td>...</td><td>...</td><td>102.3</td><td>5.7</td><td>8</td></tr><tr><td>ETAC</td><td>20.0</td><td>3.1</td><td>14</td><td>94.4</td><td>2.0</td><td>14</td></tr><tr><td>THOR</td><td>...</td><td>...</td><td>...</td><td>96.2</td><td>3.5</td><td>4</td></tr><tr><td>XFOR</td><td>18.8</td><td>1.4</td><td>5</td><td>...</td><td>...</td><td>...</td></tr><tr><td>PL8</td><td>21.9</td><td>2.0</td><td>4</td><td>...</td><td>...</td><td>...</td></tr><tr><td>ROPH</td><td>-6.3</td><td>0.3</td><td><sup>a</sup></td><td>131.0</td><td>3.0</td><td><sup>a</sup></td></tr><tr><td>CRA</td><td>-0.4</td><td>1.2</td><td>11</td><td>139.4</td><td>6.1</td><td>3</td></tr><tr><td>UCRA</td><td>-2.5</td><td>2.4</td><td>9</td><td>148.0</td><td>3.0</td><td>4</td></tr><tr><td>TAU</td><td>17.0</td><td>2.9</td><td>119</td><td>126</td><td>16</td><td>30</td></tr><tr><td>118TAU</td><td>14.7</td><td>1.1</td><td>8</td><td>112.4</td><td>5.6</td><td>6</td></tr></tbody></table>

**Note.**

<sup>a</sup> Average observables taken from Mamajek (2008).

Coma Berenices (CBER; also called Melotte 111 and Collinder 256; e.g., Casewell et al. 2006) is a massive and well-studied open cluster located at ~85 pc. Silaj & Landstreet (2014) estimate an age of $560_{-80}^{+100}$ Myr based on the Hertzsprung-Russell diagram position of its Ap-type stars. The membership lists of Casewell et al. (2006), Kraus & Hillenbrand (2007) and Casewell et al. (2014) were used here.

IC 2602 (Melotte 102; Whiteoak 1961) is one of the nearest open clusters, and is located near the Sco-Cen OB region. Its members are located at a distance of ≈150 pc (van Leeuwen 2009) and the cluster has a lithium depletion boundary age of $46_{-5}^{+6}$ Myr (Dobbie et al. 2010). The list of members published by Silaj & Landstreet (2014) and Mermilliod et al. (2009) were adopted as a starting point for BANYAN Σ.

IC 2391 (Omicron Velorum; Platais et al. 2007) is a ~50 ± 5 Myr old cluster (Barrado y Navascués et al. 2004) located at ≈150 pc, and is also in the vicinity of the Sco-Cen OB region. The membership list of Gaia Collaboration et al. (2017) was used here.

## 4.2. Associations With Partial Kinematics

This section describes the nine young associations that do not have enough known members with full 6D kinematics to build their kinematic models based on only such members. Instead, an average radial velocity or distance (or both) are adopted for the young association itself. These average distances are obtained by calculating the weighted average of all members with measured distances, where the weights are set to the inverse square of the measurement errors. The error bars on the average distance correspond to an estimate of the intrinsic distance dispersion of the members, rather than a proper measurement error of the average, and was obtained using an unbiased weighted standard deviation<fnref n="16" /> of the individual members' distance measurements, with the same weights as described above. The average radial velocities are calculated with the same method, but spectral binaries were avoided in their calculation. All average distances and/or radial velocities that were measured in this section and used in the construction of the kinematic models are listed in Table 2.

Instead of assigning the exact same average distance or radial velocity to each members with a missing observable, artificial values were drawn from a random distribution (limited to ±1σ) with a characteristic width set to the measured intrinsic dispersion of the members (listed in Table 2). This avoids artificially placing several members along lines in XYZ and UVW space—or along planes in UVW space when both radial velocity and distance are missing. The lack of full 6D kinematics for a significant number of candidates will result in a lower recovery rate of their true members by BANYAN Σ, and a larger number of contaminants from field stars. Quantifying this effect in terms of exact true-positive and false-positive rates is however not currently possible given our lack of information on the true shapes and sizes of their spatial and kinematic distributions; however Gaia-DR2 will allow us to greatly refine the kinematic models of these associations.

ε Chamaeleontis (EPSC; Mamajek et al. 2000; Feigelson et al. 2003; Murphy et al. 2013) is a young (3–5 Myr) and relatively distant (100–120 pc) association that is part of the Chamaeleon molecular cloud complex (Luhman et al. 2008). Murphy et al. (2013) refined the age of EPSC to $3.7_{-1.4}^{+4.6}$ Myr by comparing its members with the Dartmouth isochrones of Dotter et al. (2008). The membership list of Murphy et al. (2013) was adopted here.

The η Chamaeleontis cluster (ETAC; Mamajek et al. 1999) is a group of young (11 ± 3 Myr; Bell et al. 2015) stars located at a distance of ~100 pc, and in the vicinity of the Sco-Cen OB association. The membership lists of Mamajek et al. (2000) and Lyo et al. (2004) were adopted in this work.

The 32 Orionis group (THOR; Mamajek 2007; Shvonski et al. 2010; Bell et al. 2017) is a young group of ~25 Myr old stars located at ~96 pc. The age of THOR was revised to $22_{-3}^{+4}$ Myr by Bell et al. (2015) from a comparison of its members to model isochrones. Burgasser et al. (2016) recently identified the first substellar candidate member of THOR, with an estimated mass of $14 M_{Jup}$ near the planetary-mass boundary. The membership list of Bell et al. (2017) was adopted here.

The χ¹ For association (XFOR; also called Alessi 13) was identified by Dias et al. (2002), and Kharchenko et al. (2013) estimated an age of ~525 Myr based on the main-sequence turnoff. However, Mamajek (2015) argue that it could be as young as ~30 Myr due to the saturated X-ray emission of its members. Further studies will be required to address this discrepancy. Only one member of XFOR has been published with full kinematics (the triple star χ¹ For). In order to identify more members, the Dias et al. (2014) list of 4102 XFOR candidates was cross-matched with the Gaia-DR1. Of the 261 matches with a parallax measurement, only nine are located within 10 pc of the χ¹ For system in XYZ space. This indicates that the Dias et al. (2014) sample seems highly contaminated by background stars and we therefore recommend caution in its use. A literature search was performed to identify one additional radial velocity measurement for CD–37 1263, thus completing its kinematic measurements and making it a new likely bona fide member of XFOR (although its age is not investigated here). Six of the eight additional potential members are located within 5 km s⁻¹ of the star χ¹ For in UVW space if we assume the same radial velocity measurement, and they are therefore included as high-likelihood candidate members. Two additional XFOR members were identified by B. S. Alessi et al. (2018, private communication): HD 21434 and HD 17864. Both

<footnote n="16">¹⁶ See https://www.gnu.org/software/gsl/manual/html_node/Weighted-Samples.html.</footnote>
have a parallax measurement in Gaia-DR1, but only HD 17864 also has a radial velocity measurement available in the literature.

Platais 8 (or a Car; PL8) is a ~60 Myr old cluster of stars at a distance of ~130 pc identified by Platais et al. (1998). It has since received very little attention in the literature, and only four of its members benefit from full kinematic measurements (a Car, HD 76230, H Vel, OY Vel).

$\rho$ Ophiuchi (ROPH) is the nearest star-forming cloud complex to the Sun. It has been the subject of extensive studies in recent decades (e.g., see Reipurth et al. 1991; Reipurth 2008; Wilking et al. 2008). Its age is estimated at <2 Myr, and includes embedded clusters with stars believed to be as young as ~0.1 Myr (Luhman & Rieke 1999). Because this group is too distant for its members to have been directly detected by the Hipparcos mission, Mamajek (2008) used the measured parallaxes of Hipparcos stars illuminating the Lynds 1688 dark cloud, which is part of ROPH, to estimate its distance at $131 \pm 3$ pc. They also estimate an average radial velocity of $-6.3 \pm 0.3$ km s$^{-1}$ from individual radial velocity measurements of its members. The membership lists of Wilking et al. (2008) and Ducourant et al. (2017) were adopted here, and the average radial velocity and distance of Mamajek (2008) were adopted for all members with missing measurements. Only the 194 out of 340 members that have a proper motion measurement were used in the construction of the BANYAN $\Sigma$ kinematic model of ROPH. A cross-match of these with Gaia DR1 yielded 84 matches, indicating that the second data release will provide a wealth of new information on the distances and proper motions of ROPH members.

The Corona Australis (CRA) star-forming region is located at a distance of ~150 pc (Neuhäuser & Forbrich 2008; Reipurth 2008), and includes the well-studied R CrA dark cloud (e.g., see Wilking et al. 1992). Gennaro et al. (2012) estimated the age of the eclipsing binary system TY CrA between $3.8_{-0.2}^{+2.7}$ Myr and $5.2_{-0.7}^{+3.1}$ Myr, based on a comparison of the dynamical masses of its components with a set of warm and cold PISA pre-main-sequence models (Tognelli et al. 2011), respectively. Here we therefore adopt an age of ~4–5 Myr for CRA. The membership list of Neuhäuser & Forbrich (2008) was adopted here.

Several stars in the vicinity of CRA discovered by Neuhäuser et al. (2000) were found to be located between CRA and the Sco-Cen region, and at similar distances to both of these regions. This population likely constitutes of stars that formed along a filament between CRA and Sco-Cen ~10 Myr ago. We included them in the models of BANYAN $\Sigma$, and tentatively name this population Upper CrA (UCRA hereafter).

The Taurus-Auriga (TAU) star-forming region is a complex of several dark clouds located at ~130 pc, and composed of stars with ages ~1–2 Myr that share similar kinematics (e.g., see Kenyon & Hartmann 1995; Reipurth 2008). The membership lists of Luhman et al. (2009) and Esplin et al. (2014) were adopted here, without differentiating the sub-groups. Measurement errors were not provided for the TAU radial velocities measured by Wichmann et al. (2000), but they report two sets of measurements from two distinct instruments. We measured the standard deviations of the radial velocity differences for the 25 stars in their sample that were observed with both instruments, ignoring 10 spectral binaries and two stars with significantly different measurements (>8 km s$^{-1}$). We adopted

this standard deviation of 2 km s$^{-1}$ as their radial velocity measurement errors.

Mamajek (2016) identified 11 stars in the vicinity of TAU that share a larger proper motion and a closer distance to the Sun than the rest of the group (see the discussions of Currie et al. 2017 and Kraus et al. 2017). This group, named after its brightest member 118 Tau (118TAU hereafter), seems to display a slightly younger age than TAU, at ~10 Myr. The membership list of Mamajek (2016) was adopted here.

### 4.3. Rejected Associations

The Argus association was removed entirely from the models of BANYAN $\Sigma$, as Bell et al. (2015) demonstrated that it is either largely contaminated, or composed of objects that do not form a coeval association (see also Mamajek 2015). In addition, the Octans-Near (Zuckerman et al. 2013) and Hercules-Lyra associations (Gaidos 1998; Fuhrmann 2004; López-Santiago et al. 2006; Eisenbeiss et al. 2013) were not included in BANYAN $\Sigma$, as they were also demonstrated to be likely composed of non-coeval stars (Brandt et al. 2014; Mamajek 2015; Riedel et al. 2017b).

### 4.4. Discussion on Individual Objects

Individual stars that require more detailed considerations are discussed in this section. In addition to this, we note that several stars listed by different authors as bona fide members of different associations (e.g., V570 Car, CP-68 1388 and CD-69 1055) are excluded from the BANYAN $\Sigma$ kinematic models and are listed in Table 3.

AB Pic was incorrectly listed by Gagné et al. (2014) as a bona fide member of both THA and CAR. This was a consequence of Zuckerman & Song (2004) listing it as a bona fide member of THA and Torres et al. (2008) revising it to a bona fide member of CAR. Since its UVW position is at $0.7 \pm 0.7$ km s$^{-1}$ from the core of CAR and at $5.5 \pm 2.1$ km s$^{-1}$ from that of THA, here it was included in the list of CAR members (see also the discussions of Bell et al. 2015, Section B2.3, and Malo et al. 2013, Section 9.1.5).

DK Leo was listed by Malo et al. (2013) as an ambiguous candidate member between $\beta$ PMG, COL and ABDMG, because of contradictory measurements for its radial velocity (Montes et al. 2001; Kharchenko et al. 2007; López-Santiago et al. 2010), and Gagné et al. (2014) incorrectly listed it as a bona fide member of both $\beta$ PMG and CAR. Here it is excluded from the list of bona fide members until more radial velocity measurements become available.

HD 23524 was defined as a bona fide member of THA by Zuckerman et al. (2011), and Malo et al. (2013) defined it as a bona fide member of COL because it lies nearer in XYZ space, although they note that its membership is ambiguous. Here HD 23524 is excluded from the BANYAN $\Sigma$ kinematic models.

HIP 3556 was noted as a radial velocity variable and a spectroscopic double-lined binary by Kastner et al. (2017). Until a full radial velocity curve is available, this object is excluded from the BANYAN $\Sigma$ kinematic models.

2MASS J06085283-2753583 was identified by Rice et al. (2010) as a candidate member of $\beta$ PMG and it now has full kinematic measurements, but Faherty et al. (2016) showed that it is an ambiguous member, which is probably due to its small
<table><caption>**Table 3**
Objects Listed as Bona Fide Members that Were Excluded from the Construction of Kinematic Models</caption><thead><tr><th>Main Designation</th><th>Spectral Type</th><th>R.A. (hh:mm:ss)</th><th>Decl. (dd:mm:ss)</th><th>μ<sub>α</sub> cos δ (mas yr<sup>-1</sup>)</th><th>μ<sub>δ</sub> (mas yr<sup>-1</sup>)</th><th>Rad. Vel. (km s<sup>-1</sup>)</th><th>Distance (pc)</th><th>Reason for Exclusion<sup>a</sup></th><th>References</th></tr></thead><tbody><tr><td colspan="10"><strong>AB Doradus</strong></td></tr><tr><td>PX Vir</td><td>G5 V</td><td>13:03:49.46</td><td>-05:09:45.9</td><td>-191.1 ± 0.9</td><td>-218.7 ± 0.7</td><td>-9.2 ± 0.2</td><td>21.7 ± 0.4</td><td>VIS</td><td>1,2,3,2</td></tr><tr><td>2MASS J15534211-2049282</td><td>M3.4</td><td>15:53:42.09</td><td>-20:49:28.6</td><td>-10 ± 10</td><td>-20 ± 10</td><td>-7 ± 2</td><td>330 ± 80</td><td>σ<sub>kin</sub></td><td>4,5,4,4</td></tr><tr><td>HR 7214</td><td>A4 V</td><td>19:03:32.25</td><td>01:49:07.6</td><td>22.6 ± 0.2</td><td>-68.5 ± 0.3</td><td>-23 ± 2</td><td>54.9 ± 0.9</td><td>MST<sub>0</sub></td><td>6,2,6,2</td></tr><tr><td>LQ Peg</td><td>K8 V</td><td>21:31:01.86</td><td>23:20:05.2</td><td>134.6 ± 0.1</td><td>-144.9 ± 0.1</td><td>-17 ± 1</td><td>24.2 ± 0.1</td><td>VIS</td><td>1,7,8,7</td></tr><tr><td>HIP 107948</td><td>M3 Ve</td><td>21:52:10.51</td><td>05:37:33.7</td><td>106 ± 2</td><td>-147 ± 1</td><td>-15 ± 2</td><td>31 ± 5</td><td>σ<sub>kin</sub></td><td>9,9,9,2</td></tr><tr><td colspan="10"><strong>β Pictoris</strong></td></tr><tr><td>LP 353-51</td><td>M3 Ve</td><td>02:23:26.75</td><td>22:44:05.1</td><td>98.5 ± 0.2</td><td>-112.5 ± 0.1</td><td>13.0 ± 0.4</td><td>27.1 ± 0.3</td><td>VIS</td><td>1,7,10,7</td></tr><tr><td>HD 15115</td><td>F4 IV</td><td>02:26:16.33</td><td>06:17:32.4</td><td>87.97 ± 0.04</td><td>-50.35 ± 0.03</td><td>0.8 ± 0.1</td><td>48 ± 1</td><td>VIS</td><td>11,7,12,7</td></tr><tr><td>EPIC 211046195</td><td>M8.5</td><td>03:35:02.09</td><td>23:42:35.6</td><td>50 ± 10</td><td>-60 ± 10</td><td>16 ± 2</td><td>42 ± 2</td><td>VIS</td><td>4,5,4,4</td></tr><tr><td>HIP 23418 ABCD</td><td>M3 V</td><td>05:01:58.83</td><td>09:58:57.2</td><td>10 ± 10</td><td>-74 ± 6</td><td>18 ± 3</td><td>25 ± 1</td><td>VIS</td><td>13,2,14,15</td></tr></tbody></table>

**Note.**

<sup>a</sup> Reason for exclusion from the kinematic models. σ<sub>kin</sub>: Excluded from the kinematic precision constraints described in Section 5; Host: The host star was excluded; MST<sub>0</sub>: Excluded from iteration *i* of the Minimum Spanning Tree outlier rejection algorithm described in Section 5; Nσ: Excluded from the Mahalanobis distance criterion described in Section 5; VIS: Visual rejection, see Table 8.

**References.** (1) Malo et al. (2013), (2) van Leeuwen (2007), (3) Maldonado et al. (2010), (4) Shkolnik et al. (2012), (5) Monet et al. (2003), (6) Zuckerman et al. (2011), (7) Gaia Collaboration et al. (2016a), (8) Montes et al. (2001), (9) Zickgraf et al. (2005), (10) Shkolnik et al. (2017), (11) Lépine & Simon (2009), (12) Desidera et al. (2015), (13) Egret et al. (1992), (14) Song et al. (2003), (15) Riedel et al. (2014), (16) Malo et al. (2014), (17) Zacharias et al. (2010), (18) van Altena et al. (1995), (19) Houk & Swift (1999), (20) Perryman et al. (1998), (21) Lépine et al. (2013), (22) Wilson (1953), (23) White et al. (2007), (24) Griffin et al. (1988), (25) Cayrel de Strobel et al. (2001), (26) Stephenson (1986b), (27) Gray et al. (2003), (28) Nesterov et al. (1995), (29) Stephenson & Sanwal (1969), (30) Stephenson (1986a), (31) Evans (1967), (32) Royer et al. (2007), (33) Torres et al. (2006), (34) Mason et al. (2001), (35) Neuhäuser et al. (2003), (36) Zuckerman & Song (2004), (37) Gontcharov (2006), (38) Anderson & Francis (2012), (39) Bobylev et al. (2006), (40) Kraus et al. (2014), (41) Zuckerman & Webb (2000), (42) Grenier et al. (1999), (43) Pecaut & Mamajek (2013), (44) Elliott et al. (2014), (45) Webb et al. (1999), (46) Ducourant et al. (2014), (47) Pourbaix et al. (2004), (48) Kraus & Hillenbrand (2007), (49) Luhman & Steeghs (2004), (50) Shvonski et al. (2016), (51) M. J. Pecaut (2018, private communication), (52) Houk & Cowley (1975), (53) Kharchenko et al. (2007), (54) Pecaut & Mamajek (2016), (55) Houk (1978), (56) Torres et al. (2008), (57) Torres et al. (2003), (58) Zacharias et al. (2013), (59) Walter et al. (1988), (60) Kraus et al. (2017), (61) Wichmann et al. (2000), (62) Hartigan et al. (1994), (63) Welty et al. (1994), (64) Hartigan & Kenyon (2003), (65) Zacharias et al. (2012), (66) Martín et al. (2005), (67) Herbig (1977), (68) Roeser et al. (2010), (69) Hartmann et al. (1986), (70) Slesnick et al. (2006b), (71) Wichmann et al. (1997), (72) Houk (1982), (73) Erickson et al. (2011), (74) Houk & Smith-Moore (1988), (75) Cieza et al. (2007), (76) Prato (2007), (77) Guenther et al. (2007), (78) Lawrence et al. (2007), (79) Bouy & Martín (2009), (80) Donaldson et al. (2017), (81) Dahm et al. (2012), (82) Preibisch et al. (2002), (83) Mermilliod et al. (2009), (84) Kordopatis et al. (2013), (85) Galli et al. (2017), (86) Majewski et al. (2016), (87) Garcia et al. (1988). (This table is available in its entirety in machine-readable form.)

proper motion. This object was therefore not included in the kinematic models of BANYAN Σ.

## 4.5. New Bona Fide Members

New bona fide members that could be defined as such based on Gaia-DR1 data are described in this section, and are listed in Table 4.

The candidate members of THA identified by Kraus et al. (2014) were cross-matched with Gaia-DR1. Seven were found to have a parallax measurement: three of them did not match any moving group in BANYAN II and were therefore rejected (2MASS J02000918-8025009, 2MASS J02105538-4603588, and 2MASS J05332558-5117131), and the other four were confirmed as new bona fide members of THA (2MASS J02303239-4342232, 2MASS J04000382-2902165, 2MASS J04000395-2902280, 2MASS J04021648-1521297).

A similar cross-match of the Malo et al. (2014) candidate members missing a distance measurement with Gaia-DR1 yielded 21 matches, 16 of which were confirmed as new bona fide members (five in βPMG, eight in THA, two in ABDMG and one in COL). The UVW position of 2MASS J20395460 +0620118 (an ABDMG candidate from Malo et al. 2014) is a better match to the Gagné et al. (2014) position of ARG or βPMG than that of ABDMG. It was therefore categorized as

an ambiguous member until it is studied in more detail. Furthermore, Malo et al. (2014) assign 2MASS J02303239-4342232 in COL, whereas Kraus et al. (2014) call it a THA member, and it is therefore categorized as an ambiguous member in this compilation.

Shkolnik et al. (2017) performed a survey of new low-mass members in βPMG, and identified 39 new objects with signs of youth, sky position, proper motion, and radial velocities that match βPMG. As only parallaxes are still needed for them to be included in our kinematic models, their sample was cross-matched with Gaia-DR1. Five objects were found to have a parallax measurement. Four of them (HD 337919, TYC 2136-2484-1, TYC 2658-31-1, TYC 1084-672-1) have Gaia-DR1 trigonometric distances above 200 pc, preventing a credible membership in βPMG (they were also rejected as βPMG candidates by Shkolnik et al. 2017). The last object, TYC 2703-706-1, was designated as a βPMG candidate by Shkolnik et al. (2017), but has a UVW position that is located at 8.9 km s<sup>-1</sup> from the central position of βPMG, and only 3.4 km s<sup>-1</sup> of that of Columba, as defined by Gagné et al. (2014). We therefore categorize it as an ambiguous member between βPMG and COL until it is investigated further.

Six more objects in the Shkolnik et al. (2017) sample have a parallax measurement from other works in the literature
<table><caption>**Table 4**
New Bona Fide Members with Full Kinematics Compiled in this Work</caption><thead><tr><th>Designation</th><th>References<sup>a</sup></th></tr></thead><tbody><tr><td><strong>AB Doradus</strong></td><td></td></tr><tr><td>HS Psc</td><td>Malo et al. (2014)</td></tr><tr><td>HD 201919</td><td>Malo et al. (2014)</td></tr><tr><td><strong>β Pictoris</strong></td><td></td></tr><tr><td>AF Psc</td><td>Shkolnik et al. (2017)</td></tr><tr><td>2MASS J16572029-5343316</td><td>Malo et al. (2014)</td></tr><tr><td>CD-31 16041</td><td>Malo et al. (2014)</td></tr><tr><td>2MASS J19560438-3207376</td><td>Malo et al. (2014)</td></tr><tr><td>2MASS J232424896-7142211</td><td>Malo et al. (2014)</td></tr><tr><td>BD-13 6424</td><td>Malo et al. (2014)</td></tr><tr><td><strong>Columba</strong></td><td></td></tr><tr><td>GJ 1284</td><td>Malo et al. (2014)</td></tr><tr><td><strong>Tucana-Horologium</strong></td><td></td></tr><tr><td>2MASS J02303239-4342232</td><td>Kraus et al. (2014)</td></tr><tr><td>2MASS J04000382-2902165</td><td>Kraus et al. (2014)</td></tr><tr><td>2MASS J04000395-2902280</td><td>Kraus et al. (2014)</td></tr><tr><td>2MASS J04021648-1521297</td><td>Kraus et al. (2014)</td></tr><tr><td>CD-34 521</td><td>Malo et al. (2014)</td></tr><tr><td>CD-53 544</td><td>Malo et al. (2014)</td></tr><tr><td>CD-58 553</td><td>Malo et al. (2014)</td></tr><tr><td>CD-35 1167</td><td>Malo et al. (2014)</td></tr><tr><td>CD-44 1173</td><td>Malo et al. (2014)</td></tr><tr><td>2MASS J04480066-5041255</td><td>Malo et al. (2014)</td></tr><tr><td>2MASS J05332558-5117131</td><td>Malo et al. (2014)</td></tr><tr><td>2MASS J23261069-7323498</td><td>Malo et al. (2014)</td></tr><tr><td><strong>χ<sup>1</sup> For</strong></td><td></td></tr><tr><td>CD-37 1263<sup>b</sup></td><td>Dias et al. (2014)</td></tr><tr><td>HD 17864</td><td>B. S. Alessi (2018, private communication)</td></tr></tbody></table>

**Notes.** This table lists objects that were designated as candidate members and that we confirm as bona fide members with full kinematics from compiling their missing measurements. See Section 4 for more details.

<sup>a</sup> Reference that designated this object as a candidate member of the young association.

<sup>b</sup> The age of CD-37 1263 has not been investigated here; verifying its youth is still necessary to confirm that it is a true bona fide member.

(Shkolnik et al. 2012; Riedel et al. 2014); five/six were already included in the list of bona fide members presented here. The remaining star, AF Psc, has a parallax measurement by van Altena et al. (1995), which seems to have been overlooked in previous studies. It was thus added to the list of bona fide members of βPMG.

Riedel et al. (2017a) presented several new M-type young moving group candidates; however none of them benefit from a parallax distance measurement either in Gaia-DR1 or elsewhere in the literature, hence they were not included in the list of bona fide members.

## 4.6. Calculation of the 6D Kinematics

The Galactic positions XYZ and space velocities UVW, in a right-handed system where U points toward the Galactic center, were calculated for all members by assuming Gaussian error bars in sky position, proper motion, radial velocity and parallax. A 10⁴-element Monte Carlo approach was used to propagate error bars in XYZUVW space, by adopting the standard deviation of each coordinate as its measurement error, and therefore assuming that error bars are Gaussian in XYZUVW space. The resulting list of bona fide members is presented in Table 5 with the corresponding references in

Table 6, and their various designations are listed in Table 7. A list of new bona fide members confirmed in this work is given in Table 4, and their positions on the sky are displayed in Figure 1.

## 5. Kinematic Models of Young Associations

The bona fide members compiled in Table 5 were used to build the kinematic models of young associations considered in BANYAN Σ. Objects with total error bars on their Galactic position (√σ<sup>2</sup><sub>X</sub> + σ<sup>2</sup><sub>Y</sub> + σ<sup>2</sup><sub>Z</sub>) above 20% of their distance or on their space velocity (√σ<sup>2</sup><sub>U</sub> + σ<sup>2</sup><sub>V</sub> + σ<sup>2</sup><sub>W</sub>) above 8 km s<sup>-1</sup> were excluded. Companions in binary systems were ignored, and the spatial-kinematic position of binary systems was approximated as that of the primary star to avoid the need for model-dependent mass estimates.

A rejection algorithm based on minimum spanning trees (MSTs; e.g., see Allison et al. 2009; Gagné et al. 2015c) was used to ignore outliers in spatial and kinematic space. Groups that required adopting average radial velocities or distances for some high-likelihood members were exempted from this rejection step because of their small number of bona fide members (118TAU, CRA, EPSC, ETAC, PL8, ROPH, TAU, THOR, UCRA, and XFOR). If these groups are contaminated by outliers that originate from larger spatial or kinematic distributions, this exemption will result in models that are artificially biased to larger sizes, therefore increasing the rate of contamination that these groups may be subject to. The discovery of more members will be necessary to better assess and correct this effect. A spanning tree is built by connecting each star of a young association in XYZ or UVW space with straight lines while avoiding loops; the MST is the spanning tree with the shortest total length. MSTs provide a measurement of scale that does not depend on the shape of a distribution, and therefore does not require making the assumption that the stars are normally distributed, or aligned with the XYZUVW axes.

For each association containing a total of N<sub>k</sub> members, the algorithm of Cartwright & Whitworth (2004) was used to build N<sub>k</sub>+1 MSTs, in both XYZ and UVW spaces separately. The first spatial and kinematic MSTs with respective total lengths L<sub>spa</sub> and L<sub>kin</sub> include all members, and the N<sub>k</sub> additional MSTs ignore one member at a time, and have respective total lengths L<sub>i,spa</sub> and L<sub>i,kin</sub>.

Relative lengths L<sub>i,rel</sub> of MSTs ignoring one star each were then calculated, and an arbitrary threshold L<sub>t</sub> was set at a value 5% smaller than the 90% percentile value ⟨L<sub>i,rel</sub>⟩<sub>90%</sub> of all the relative MST lengths:

$$L_{i,rel} = \sqrt{\left(\frac{L_{i,spa}}{L_{spa}}\right)^2 + \left(\frac{L_{i,kin}}{L_{kin}}\right)^2}, \\ L_t = 0.95 \langle L_{i,rel} \rangle_{90\%}.$$

Any individual MST with relative length L<sub>i,rel</sub> < L<sub>t</sub> indicates that removing a single star i significantly shortens the spatial and/or kinematic size of the young association, and therefore that star i is an outlier. Such outliers were removed, and an iterative rejection process with a more conservative arbitrary threshold L<sub>t</sub> = 0.9⟨L<sub>i,rel</sub>⟩<sub>90%</sub> was subsequently used until no stars are rejected.

In a few cases, the bona fide members that survived the MST rejection displayed a clumped distribution of members with
<table><caption>**Table 5**
Literature Compilation of Bona Fide Members</caption><thead><tr><th>Main<br>Designation</th><th>Spectral<br>Type</th><th>R.A.<br>(hh:mm:ss)</th><th>Decl.<br>(dd:mm:ss)</th><th>μ<sub>α</sub> cos δ<br>(mas yr<sup>-1</sup>)</th><th>μ<sub>δ</sub><br>(mas yr<sup>-1</sup>)</th><th>Rad. Vel.<br>(km s<sup>-1</sup>)</th><th>Distance<br>(pc)</th><th>References</th></tr></thead><tbody><tr><td colspan="9"><em>AB Doradus</em></td></tr><tr><td>2MASS J00192626+4614078</td><td>M8 β</td><td>00:19:26.26</td><td>46:14:07.8</td><td>119.4 ± 0.9</td><td>-75.4 ± 0.9</td><td>-20 ± 3</td><td>39 ± 2</td><td>1, 2, 3, 2</td></tr><tr><td>BD+54 144 A</td><td>F8 V</td><td>00:45:51.06</td><td>54:58:39.1</td><td>96.40 ± 0.03</td><td>-73.97 ± 0.04</td><td>-15 ± 2</td><td>50.3 ± 0.9</td><td>4, 5, 6, 5</td></tr><tr><td>— BD+54 144 B</td><td>K3</td><td>00:45:51.23</td><td>54:58:40.8</td><td>...</td><td>...</td><td>...</td><td>...</td><td>7, ..., ..., ...</td></tr><tr><td>2MASS J00470038+6803543</td><td>L6–L8 γ</td><td>00:47:00.39</td><td>68:03:54.4</td><td>385 ± 1</td><td>-201 ± 1</td><td>-20 ± 1</td><td>12.2 ± 0.3</td><td>..., 2, 8, 2</td></tr><tr><td>G 132–51 B</td><td>M2.6</td><td>01:03:42.23</td><td>40:51:13.6</td><td>132 ± 5</td><td>-164 ± 5</td><td>-10.6 ± 0.3</td><td>30 ± 2</td><td>9, 9, 9, 9</td></tr><tr><td>— G 132–50</td><td>M0</td><td>01:03:40.30</td><td>40:51:26.7</td><td>126.9 ± 0.1</td><td>-166.18 ± 0.09</td><td>...</td><td>33.0 ± 0.6</td><td>10, 5, ..., 5</td></tr><tr><td>— G 132–51 C</td><td>M3.8</td><td>01:03:42.44</td><td>40:51:13.3</td><td>...</td><td>...</td><td>-10.9 ± 0.4</td><td>...</td><td>9, ..., 9, ...</td></tr><tr><td>HIP 6276</td><td>G0 V</td><td>01:20:32.38</td><td>-11:28:05.8</td><td>111.44 ± 0.06</td><td>-136.95 ± 0.05</td><td>8.3 ± 0.4</td><td>35.1 ± 0.4</td><td>11, 5, 12, 5</td></tr><tr><td>G 269–153 A</td><td>M4.3</td><td>01:24:27.85</td><td>-33:55:11.4</td><td>180 ± 20</td><td>-110 ± 20</td><td>19 ± 3</td><td>25 ± 1</td><td>9, 13, 9, 9</td></tr><tr><td>— G 269–153 B</td><td>M4.6</td><td>01:24:27.96</td><td>-33:55:09.9</td><td>...</td><td>...</td><td>18 ± 1</td><td>...</td><td>9, ..., 9, ...</td></tr></tbody></table>

**Note.** References—The references to this table are listed in Table 6.

(This table is available in its entirety in machine-readable form.)
<table><caption>**Table 6**
References for Literature Compilation of Bona Fide Members</caption><thead><tr><th>References</th></tr></thead><tbody><tr><td>(1) Gagné et al. (2015b), (2) Liu et al. (2016), (3) Reiners & Basri (2009), (4) Jaschek et al. (1964), (5) Gaia Collaboration et al. (2016a), (6) Bobylev et al. (2006), (7) Zuckerman & Song (2004), (8) Faherty et al. (2016), (9) Shkolnik et al. (2012), (10) Lépine et al. (2013), (11) Malo et al. (2013), (12) White et al. (2007), (13) Monet et al. (2003), (14) Malo et al. (2014), (15) Schlieder et al. (2010), (16) van Leeuwen (2007), (17) Ivanov (2008), (18) Tokovinin & Smekhov (2002), (19) Garrison & Gray (1994), (20) Gontcharov (2006), (21) Bystrov et al. (1994), (22) Egret et al. (1992), (23) Anderson & Francis (2012), (24) Houk & Cowley (1975), (25) Zuckerman et al. (2011), (26) Blake et al. (2010), (27) Chubak & Marcy (2011), (28) Kordopatis et al. (2013), (29) Høg et al. (2000), (30) Torres et al. (2006), (31) Zacharias et al. (2004a), (32) Gray et al. (2006), (33) Martín & Brandner (1995), (34) Kharchenko et al. (2007), (35) Messina et al. (2010), (36) Montes et al. (2001), (37) Reid et al. (2004), (38) Kunder et al. (2017), (39) Terrien et al. (2015), (40) Knapp et al. (2004), (41) Dupuy & Liu (2012), (42) Gagné et al. (2015a), (43) Dieterich et al. (2014), (44) Röser et al. (2008), (45) Abt & Morrell (1995), (46) Zuckerman et al. (2001a), (47) Zacharias et al. (2010), (48) Riedel et al. (2014), (49) Valenti & Fischer (2005), (50) Holmberg et al. (2007), (51) Song et al. (2003), (52) Macintosh et al. (2015), (53) Zacharias et al. (2004b), (54) Allers & Liu (2013), (55) Lépine & Simon (2009), (56) Bobylev & Bajkova (2007), (57) Gizis et al. (2002), (58) Bonnefoy et al. (2013), (59) Torres et al. (2009), (60) Kiss et al. (2011), (61) Corbally (1984), (62) Liu et al. (2013), (63) Allers et al. (2016), (64) Shkolnik et al. (2017), (65) Gray et al. (2003), (66) Eggl et al. (2013), (67) King et al. (2003), (68) Levato & Abt (1978), (69) Evans (1967), (70) Fabricius et al. (2002), (71) Gray & Garrison (1987), (72) Mamajek et al. (2010), (73) Zuckerman et al. (2006), (74) Salim & Gould (2003), (75) Nordström et al. (2004), (76) Desidera et al. (2015), (77) Koen et al. (2010), (78) Griffin et al. (1988), (79) Joy & Wilson (1949), (80) Wilson (1953), (81) Perryman et al. (1998), (82) Stephenson (1986a), (83) Hussain et al. (2006), (84) Nesterov et al. (1995), (85) Cayrel de Strobel et al. (2001), (86) Stephenson (1986b), (87) Cenarro et al. (2009), (88) Gebran et al. (2010), (89) Stocke et al. (1991), (90) Gray et al. (2001), (91) Karataş et al. (2004), (92) Cenarro et al. (2007), (93) Pourbaix et al. (2004), (94) Morgan & Hiltner (1965), (95) Keenan & McNeil (1989), (96) Bilfková et al. (2010), (97) Morgan & Keenan (1973), (98) Paunzen et al. (2001), (99) Benedict et al. (2014), (100) Mermilliod et al. (2009), (101) Cowley & Fraquelli (1974), (102) Morse et al. (1991), (103) de Bruijne & Eilers (2012), (104) Christy & Walker (1969), (105) van Belle & von Braun (2009), (106) Tomkin et al. (1995), (107) Adams et al. (1935), (108) Gray & Garrison (1989a), (109) Kraft (1965), (110) Patel et al. (2013), (111) Abt & Levy (1985), (112) Wilson (1962), (113) Maderak et al. (2013), (114) Royer et al. (2007), (115) Nassau & Macrae (1955), (116) Cowley et al. (1969), (117) Zuckerman & Webb (2000), (118) Cruz et al. (2007), (119) Levato (1975), (120) Gagné et al. (2015c), (121) Głěbocki & Gnaciński (2005), (122) Kraus et al. (2014), (123) Moór et al. (2006), (124) Houk (1982), (125) Ducourant et al. (2014), (126) Elliott et al. (2014), (127) Gagné et al. (2017b), (128) Teixeira et al. (2008), (129) Pecaut & Mamajek (2013), (130) Weinberger et al. (2013), (131) Webb et al. (1999), (132) Torres et al. (2003), (133) Shkolnik et al. (2011), (134) Looper et al. (2010b), (135) Donaldson et al. (2016), (136) Looper et al. (2010a), (137) Schneider et al. (2012), (138) Mohanty et al. (2003), (139) Zacharias et al. (2013), (140) Gizis et al. (2007), (141) Rodriguez et al. (2011), (142) Mamajek (2005), (143) Kraus & Hillenbrand (2007), (144) Yoss & Griffin (1997), (145) Massarotti et al. (2008), (146) Lyo et al. (2004), (147) Girard et al. (2011), (148) Luhman & Steeghs (2004), (149) Lopez Martí et al. (2013), (150) Bell et al. (2017), (151) Zacharias et al. (2015), (152) Shvonski et al. (2016), (153) Roeser et al. (2010), (154) Alcalá et al. (2000), (155) Mace et al. (2009), (156) Edwards (1976), (157) Houk (1978), (158) Jackson & Stoy (1955), (159) Mamajek (2015), (160) Murphy et al. (2013), (161) Torres et al. (2008), (162) Terranegra et al. (1999), (163) Kastner et al. (2012), (164) Guenther et al. (2007), (165) Mamajek et al. (2000), (166) Grenier et al. (1999), (167) Hales et al. (2014), (168) Grady et al. (2004), (169) Luhman (2004), (170)</td></tr></tbody></table>

<table><caption>**Table 6**
(Continued)</caption><thead><tr><th>References</th></tr></thead><tbody><tr><td>Riaz et al. (2006), (171) E. Bubar et al. (2018, in preparation), (172) Covino et al. (1997), (173) Li & Hu (1998), (174) Kraus et al. (2017), (175) Abt (2008), (176) Abt (2004), (177) Slesnick et al. (2006b), (178) Hiltner et al. (1969), (179) Pecaut & Mamajek (2016), (180) Song et al. (2012), (181) Chen et al. (2011), (182) Levenhagen & Leister (2006), (183) Lutz & Lutz (1977), (184) Moór et al. (2011), (185) Walter et al. (1988), (186) Riviere-Marichalar et al. (2012), (187) Wichmann et al. (2000), (188) Herbig (1977), (189) White & Basri (2003), (190) Hartmann et al. (1987), (191) Sartoretti et al. (1998), (192) Alves de Oliveira et al. (2012), (193) Xiao et al. (2012), (194) Hartmann et al. (1986), (195) Patterer et al. (1993), (196) Donati et al. (1997), (197) Herbig et al. (1986), (198) Hessman & Guenther (1997), (199) Appenzeller et al. (1988), (200) Mundt et al. (1983), (201) Sestito et al. (2008), (202) Hartigan & Kenyon (2003), (203) Strassmeier (2009), (204) Hartigan et al. (1994), (205) Monin et al. (2010), (206) Herbig (1990), (207) Mathieu et al. (1997), (208) Joy (1949), (209) Zacharias et al. (2012), (210) Mooley et al. (2013), (211) Esplin et al. (2014), (212) Cohen & Kuhi (1979), (213) Mora et al. (2001), (214) Duchêne et al. (1999), (215) Gahm et al. (1999), (216) Bragança et al. (2012), (217) Hube (1970), (218) Buscombe (1969), (219) Krautter et al. (1997), (220) Mamajek et al. (2002), (221) Köhler et al. (2000), (222) Erickson et al. (2011), (223) Cieza et al. (2007), (224) Prato (2007), (225) Struve & Rudkjøbing (1949), (226) Bouvier & Appenzeller (1992), (227) Wilking et al. (2005), (228) Lawrence et al. (2007), (229) Ansdell et al. (2016), (230) Manara et al. (2015), (231) Kurosawa et al. (2006), (232) Ricci et al. (2010), (233) Ducourant et al. (2017), (234) Suárez et al. (2006), (235) Martín et al. (1998), (236) Slesnick et al. (2006a), (237) Lodieu (2013), (238) Rydgren (1980), (239) Orellana et al. (2012), (240) Houk & Smith-Moore (1988), (241) Dahm et al. (2012), (242) Preibisch et al. (1998), (243) Donaldson et al. (2017), (244) Siebert et al. (2011), (245) Bouy & Martín (2009), (246) Preibisch & Zinnecker (1999), (247) Walter et al. (1994), (248) Abt (1981), (249) Preibisch et al. (2001), (250) Cucchiaro et al. (1976), (251) Carpenter et al. (2006), (252) Beavers & Cook (1980), (253) Donati et al. (2006), (254) Abt (2009), (255) Galli et al. (2017), (256) Majewski et al. (2016), (257) Cottaar et al. (2015), (258) McCarthy & Treanor (1964), (259) Breger (1984), (260) Abt & Levato (1978), (261) Mendoza V (1956), (262) Gray & Garrison (1989b), (263) Binnendijk (1946), (264) Walter et al. (1997), (265) Zacharias et al. (2017), (266) Melo (2003), (267) Carmona et al. (2007), (268) Forbrich & Preibisch (2007), (269) Vieira et al. (2003), (270) Corporon et al. (1996), (271) Garcia et al. (1988), (272) Messina et al. (2011), (273) De Silva et al. (2013), (274) Pickles & Depagne (2010), (275) Moór et al. (2013), (276) Bourgés et al. (2014), (277) Anderson & Francis (2012), (278) Bobylev et al. (2006), (279) Kharchenko et al. (2007), (280) White et al. (2007), (281) Torres et al. (2006), (282) Siebert et al. (2011), (283) Chubak et al. (2012), (284) Mermilliod et al. (2009).</td></tr></tbody></table>

some remaining outliers that were visually easy to recognize. They typically survived the MST selection cuts because they are located near at least one other outlier. As a consequence, the 1D projections of the multivariate Gaussian models compared to the location of bona fide members (see Figure 2) were visually inspected to impose additional rejection criteria on the distribution of true members. These criteria are listed in Table 8.

The bona fide members that were rejected from the model construction are listed in Table 3. An average of three objects (typically less than eight objects) were rejected in each young association. Only the PLE and HYA had more rejected members (10 and 15, respectively), but given their large number of members this represents less than 8% of their
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Sky distribution of young association members that were used here to build the models of BANYAN Σ. The Galactic plane ($|b| < 15°$) is designated with the gray region. The nearest young associations cover much larger fractions of the sky, which makes it harder to recognize their members without measuring their full 6D kinematics. Most of the young association members are located in the Southern hemisphere, with a few notable exceptions (UMA, CBER, PLE, HYA, TAU, and 118TAU). See Section 4 for more details.</figcaption></figure>

<table><caption>**Table 7**
List of Bona Fide Members' Designations</caption><thead><tr><th>Main</th><th>2MASS</th><th>AllWISE</th><th>Gaia</th><th>Other</th></tr></thead><tbody><tr><td><em>AB Doradus</em></td><td></td><td></td><td></td><td></td></tr><tr><td>2MASS J00192626+4614078</td><td>J00192626+4614078</td><td>J001926.39+461406.8</td><td>...</td><td>...</td></tr><tr><td>BD+54 144 A</td><td>J00455088+5458402</td><td>J004551.02+545839.6</td><td>417565757132068096</td><td>HD 4277 A, HIP 3589 A</td></tr><tr><td>— BD+54 144 B</td><td>...</td><td>...</td><td>...</td><td>HD 4277 B, HIP 3589 B</td></tr><tr><td>2MASS J00470038+6803543</td><td>J00470038+6803543</td><td>J004701.09+680352.2</td><td>529737830321549952</td><td>...</td></tr><tr><td>G 132–51 B</td><td>J01034210+4051158</td><td>J010342.24+405114.2</td><td>374400889126932096</td><td>...</td></tr><tr><td>— G 132–50</td><td>J01034013+4051288</td><td>J010340.24+405127.4</td><td>374400957846408192</td><td>...</td></tr><tr><td>— G 132–51 C</td><td>...</td><td>...</td><td>374400893422315648</td><td>...</td></tr><tr><td>HIP 6276</td><td>J01203226–1128035</td><td>J012032.34–112805.2</td><td>2470272808484339200</td><td>CD–12 243</td></tr><tr><td>G 269–153 A</td><td>J01242767–3355086</td><td>J012427.84–335510.0</td><td>5015892473055253120</td><td>...</td></tr><tr><td>— G 269–153 B</td><td>...</td><td>...</td><td>5015892473055253248</td><td>...</td></tr></tbody></table>

(This table is available in its entirety in machine-readable form.)

<table><caption>**Table 8**
Visual Selection Cuts Applied to Moving Groups with Noticeable Outliers</caption><thead><tr><th>Association</th><th>Rejection Criterion</th></tr></thead><tbody><tr><td>ABDMG</td><td>V &gt; -24 km s<sup>-1</sup></td></tr><tr><td>βPMG</td><td>V &gt; -13 km s<sup>-1</sup><br>W &gt; -4 km s<sup>-1</sup></td></tr><tr><td>CBER</td><td>W &lt; -3 km s<sup>-1</sup></td></tr><tr><td>ETAC</td><td>Z &lt; -37 pc</td></tr><tr><td>HYA</td><td>V &lt; -22 km s<sup>-1</sup></td></tr><tr><td>LCC</td><td>U &lt; -15 km s<sup>-1</sup></td></tr><tr><td>PLE</td><td>U &gt; 0 km s<sup>-1</sup></td></tr><tr><td>THA</td><td>Z &gt; -20 pc<br>U &lt; 12.5 km s<sup>-1</sup><br>-4 km s<sup>-1</sup> &lt; W &lt; 2 km s<sup>-1</sup></td></tr><tr><td>THOR</td><td>U &lt; -25 km s<sup>-1</sup></td></tr><tr><td>USCO</td><td>U &gt; 10 km s<sup>-1</sup></td></tr></tbody></table>

populations. These objects could still be true members of their respective associations with low-quality or inaccurate kinematic measurements, therefore we do not consider that they are necessarily non-members. The exclusion of such outliers will avoid biasing the spatial and kinematic sizes of the BANYAN Σ models to artificially large values, which would result in larger rates of contamination, as long as they are either non-members or suffer from inaccurate or low-quality kinematic measurements. If some of them are true members of a young association, it is likely that other objects with similar kinematics exist and are not currently known. In such a case, the gradual discovery of moving group members at the edges of the current models with BANYAN Σ, or other methods to identify new moving groups (e.g., see Oh et al. 2017) will make it possible to uncover them. Careful searches using BANYAN Σ with the kinematics-only mode described in
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Multivariate Gaussian model of TWA. The 1, 2 and 3σ projected contours of the multivariate Gaussian model are displayed as orange lines, and black points represent individual bona fide members. The residuals resulting from the difference of a 2D kernel density estimate distribution using Silverman's rule of thumb and the multivariate Gaussian models are displayed in the background of each 2D projection. Green shades indicate an over-density of members compared to the model, and blue shades indicate an under-density. Unidimensional distributions (i.e., histograms) of the bona fide members are displayed as green bars. The thick black line represents a 1D kernel density estimate using Silverman's rule, and the thick orange line represents the projection of the multivariate Gaussian model. In the tridimensional projection figures, a single 1σ contour of the multivariate Gaussian model is displayed. Projections of the bona fide members' positions on the three axis planes are displayed with blue spheres to facilitate viewing. The complete figure set (27 images), one for each young association, is available in the online journal. See Section 3.1 for more details.</figcaption></figure>

(The complete figure set (27 images) is available.)

Section 3.9 will also make it possible to identify such groups of previously unrecognized members that are outside of the spatial dimensions of the current models, but not outside of their kinematic dimensions. Such investigations are left for

future work, and all members that are rejected here will be ignored in what follows.

The weighted averages $Q_{,\text{avg}}$ and covariances $C_{ij}$ of all spatial-kinematic coordinates $\{Q_i\}$ were calculated, where the
<table><caption>**Table 9**

Parameters for the Central Location and Variances of the Multivariate Gaussian Models of Young Associations</caption><thead><tr><th rowspan="2">Asso.</th><th>⟨X⟩</th><th>⟨Y⟩</th><th>⟨Z⟩</th><th>⟨U⟩</th><th>⟨V⟩</th><th>⟨W⟩</th><th>Σ<sup>1/2</sup><sub>00</sub></th><th>Σ<sup>1/2</sup><sub>11</sub></th><th>Σ<sup>1/2</sup><sub>22</sub></th><th>Σ<sup>1/2</sup><sub>33</sub></th><th>Σ<sup>1/2</sup><sub>44</sub></th><th>Σ<sup>1/2</sup><sub>55</sub></th></tr><tr><th colspan="2">(pc)</th><th colspan="2">(km s<sup>-1</sup>)</th><th colspan="2">(pc)</th><th colspan="5">(km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>118TAU</td><td>-102.3</td><td>-4.8</td><td>-9.9</td><td>-12.8</td><td>-19.1</td><td>-9.2</td><td>12.7</td><td>2.4</td><td>1.8</td><td>2.1</td><td>2.8</td><td>1.6</td></tr><tr><td>ABDMG</td><td>-6.0</td><td>-7.2</td><td>-8.8</td><td>-7.2</td><td>-27.6</td><td>-14.2</td><td>21.4</td><td>20.3</td><td>16.3</td><td>1.4</td><td>1.0</td><td>1.8</td></tr><tr><td>βPMG</td><td>4.1</td><td>-6.7</td><td>-15.7</td><td>-10.9</td><td>-16.0</td><td>-9.0</td><td>29.3</td><td>14.0</td><td>9.0</td><td>2.2</td><td>1.2</td><td>1.0</td></tr><tr><td>CAR</td><td>6.7</td><td>-50.5</td><td>-15.5</td><td>-10.66</td><td>-21.92</td><td>-5.48</td><td>10.0</td><td>18.1</td><td>12.6</td><td>0.67</td><td>1.02</td><td>1.01</td></tr><tr><td>CARN</td><td>0.7</td><td>-28.1</td><td>-4.3</td><td>-25.3</td><td>-18.1</td><td>-2.3</td><td>7.8</td><td>20.8</td><td>17.3</td><td>3.2</td><td>1.9</td><td>2.0</td></tr><tr><td>CBER</td><td>-6.0</td><td>-5.1</td><td>84.9</td><td>-2.30</td><td>-5.51</td><td>-0.61</td><td>3.3</td><td>3.3</td><td>4.5</td><td>0.53</td><td>0.44</td><td>0.71</td></tr><tr><td>COL</td><td>-25.9</td><td>-25.9</td><td>-21.4</td><td>-11.90</td><td>-21.28</td><td>-5.66</td><td>12.1</td><td>23.0</td><td>17.8</td><td>1.04</td><td>1.29</td><td>0.75</td></tr><tr><td>CRA</td><td>132.45</td><td>-0.21</td><td>-42.43</td><td>-3.7</td><td>-15.7</td><td>-8.8</td><td>3.71</td><td>0.75</td><td>2.04</td><td>1.3</td><td>2.2</td><td>2.2</td></tr><tr><td>EPSC</td><td>49.9</td><td>-84.8</td><td>-25.6</td><td>-9.9</td><td>-19.3</td><td>-9.7</td><td>2.5</td><td>3.6</td><td>4.0</td><td>1.6</td><td>2.2</td><td>2.0</td></tr><tr><td>ETAC</td><td>33.65</td><td>-81.36</td><td>-34.81</td><td>-10.0</td><td>-22.3</td><td>-11.7</td><td>0.65</td><td>0.98</td><td>0.71</td><td>1.6</td><td>2.8</td><td>1.8</td></tr><tr><td>HYA</td><td>-38.5</td><td>0.8</td><td>-15.8</td><td>-42.27</td><td>-18.79</td><td>-1.47</td><td>7.4</td><td>4.4</td><td>2.9</td><td>2.01</td><td>0.94</td><td>1.10</td></tr><tr><td>IC2391</td><td>1.9</td><td>-148.1</td><td>-18.0</td><td>-23.04</td><td>-14.89</td><td>-5.48</td><td>1.3</td><td>6.4</td><td>1.4</td><td>1.10</td><td>3.40</td><td>0.78</td></tr><tr><td>IC2602</td><td>47.4</td><td>-137.6</td><td>-12.6</td><td>-8.22</td><td>-20.60</td><td>-0.58</td><td>1.5</td><td>5.4</td><td>1.1</td><td>1.18</td><td>2.61</td><td>0.65</td></tr><tr><td>LCC</td><td>54.3</td><td>-94.2</td><td>5.8</td><td>-7.8</td><td>-21.5</td><td>-6.2</td><td>11.9</td><td>12.4</td><td>13.7</td><td>2.7</td><td>3.8</td><td>1.8</td></tr><tr><td>OCT</td><td>4.0</td><td>-96.9</td><td>-59.7</td><td>-13.7</td><td>-3.3</td><td>-10.1</td><td>78.3</td><td>25.8</td><td>8.8</td><td>2.4</td><td>1.3</td><td>1.4</td></tr><tr><td>PL8</td><td>10.6</td><td>-124.5</td><td>-13.9</td><td>-11.01</td><td>-22.89</td><td>-3.59</td><td>7.0</td><td>11.6</td><td>4.5</td><td>1.15</td><td>1.96</td><td>0.74</td></tr><tr><td>PLE</td><td>-118.9</td><td>28.5</td><td>-54.4</td><td>-6.7</td><td>-28.0</td><td>-14.0</td><td>7.7</td><td>3.5</td><td>4.2</td><td>1.7</td><td>1.8</td><td>1.2</td></tr><tr><td>ROPH</td><td>124.79</td><td>-15.23</td><td>37.60</td><td>-5.9</td><td>-13.5</td><td>-7.9</td><td>1.33</td><td>0.51</td><td>0.66</td><td>1.3</td><td>4.7</td><td>4.3</td></tr><tr><td>TAU</td><td>-116.3</td><td>6.7</td><td>-35.9</td><td>-14.3</td><td>-9.3</td><td>-8.8</td><td>11.4</td><td>10.8</td><td>10.1</td><td>3.1</td><td>4.5</td><td>3.4</td></tr><tr><td>THA</td><td>5.4</td><td>-20.1</td><td>-36.1</td><td>-9.79</td><td>-20.94</td><td>-0.99</td><td>19.4</td><td>12.4</td><td>3.8</td><td>0.87</td><td>0.79</td><td>0.72</td></tr><tr><td>THOR</td><td>-88.4</td><td>-25.7</td><td>-23.9</td><td>-12.8</td><td>-18.8</td><td>-9.0</td><td>4.1</td><td>6.9</td><td>5.1</td><td>2.2</td><td>2.2</td><td>2.0</td></tr><tr><td>TWA</td><td>14.4</td><td>-47.7</td><td>22.7</td><td>-11.6</td><td>-17.9</td><td>-5.6</td><td>12.2</td><td>9.7</td><td>3.9</td><td>1.8</td><td>1.8</td><td>1.6</td></tr><tr><td>UCL</td><td>107.5</td><td>-60.9</td><td>26.5</td><td>-4.7</td><td>-19.7</td><td>-5.2</td><td>21.0</td><td>19.6</td><td>13.5</td><td>3.8</td><td>3.0</td><td>1.7</td></tr><tr><td>UCRA</td><td>142.1</td><td>-1.2</td><td>-39.2</td><td>-3.7</td><td>-17.1</td><td>-8.0</td><td>7.3</td><td>2.4</td><td>5.9</td><td>3.0</td><td>1.8</td><td>1.2</td></tr><tr><td>UMA</td><td>-7.5</td><td>9.9</td><td>21.9</td><td>14.8</td><td>1.8</td><td>-10.2</td><td>3.1</td><td>1.5</td><td>1.1</td><td>1.0</td><td>1.2</td><td>2.6</td></tr><tr><td>USCO</td><td>121.2</td><td>-17.0</td><td>48.9</td><td>-4.9</td><td>-14.2</td><td>-6.5</td><td>17.0</td><td>8.2</td><td>8.9</td><td>3.7</td><td>3.2</td><td>2.3</td></tr><tr><td>XFOR</td><td>-27.1</td><td>-46.3</td><td>-84.2</td><td>-12.54</td><td>-22.24</td><td>-6.26</td><td>4.7</td><td>3.8</td><td>4.4</td><td>0.96</td><td>1.41</td><td>2.21</td></tr></tbody></table>

Note. The average Galactic positions ⟨X⟩, ⟨Y⟩ and ⟨Z⟩ and space velocities ⟨U⟩, ⟨V⟩ and ⟨W⟩ correspond to the components of the $\bar{\tau}$ vector defining the center of the multivariate Gaussian model. Parameters $\Sigma_{00}$ through $\Sigma_{55}$ represent the diagonal elements of the covariance matrix, and correspond to the dispersion of Galactic positions and space velocities along the principal axes of the multivariate Gaussian, which are not necessarily aligned with the Galactic coordinates axes. See Section 5 for more details. The covariances of the multivariate Gaussian models are given in a FITS file format with the BANYAN $\Sigma$ codes (Gagné, 2018a, 2018b).

weight of a given star is set to the inverse square of its spatial and kinematic error bars $\sigma_{k,spa}$ and $\sigma_{k,kin}$ relative to the association averages $\sigma_{avg,spa}$ and $\sigma_{avg,kin}$, added in quadrature:

$$Q_{i,\text{avg}} = \frac{\sum_k w_k Q_i}{w_{\text{tot}}},$$ $$C_{ij} = \frac{\sum_k w_k (Q_i - Q_{i,\text{avg}})(Q_j - Q_{j,\text{avg}})}{w_{\text{tot}}(w_{\text{tot}}^2 - \sum_k w_k^2)(\sum_k w_k^2)^{-1}},$$ $$w_k = \left( \left( \frac{\sigma_{k,\text{spa}}}{\sigma_{\text{avg,spa}}} \right)^2 + \left( \frac{\sigma_{k,\text{kin}}}{\sigma_{\text{avg,kin}}} \right)^2 \right)^{-1},$$ $$w_{\text{tot}} = \sum_k w_k.$$

The values of weights were set to a maximum of $w_k < 50$, corresponding to measurement errors 10 times more precise than the association average, to avoid the possibility of a very small number of precise XYZUVW measurements bearing too much weight in the kinematic models.

The covariance matrix $\bar{\Sigma}$ and center vector $\bar{\tau}$ of an association were then built from its components $Q_{i,avg}$ and $C_{ij}$, and the covariance matrix was regularized<fnref n="17" /> to avoid

numerical problems. This was done through a singular value decomposition of the covariance matrix:

$$\bar{\Sigma} = \bar{\Sigma}_U \bar{\Sigma}_{sv} \bar{\Sigma}_V^T,$$

where $\bar{\Sigma}_{sv}$ is a diagonal matrix containing the singular values. If the determinants $|\bar{\Sigma}_U|$ or $|\bar{\Sigma}_V|$ were found to be negative, random noise with a standard deviation equal to half the error bars was added to the XYZUVW coordinates of all members until $\bar{\Sigma}$ was found to be nonsingular. If an association contained less than 30 members, the three spatial and three kinematic singular values were forced to a minimum of 1 pc and 0.2 km s⁻¹, respectively.

Because the addition of noise in the regularization process can break the symmetry of the covariance matrix, the non-diagonal elements of the covariance matrix are forced to be symmetric by setting the values of both $\Sigma_{ij}$ and $\Sigma_{ji}$ to the average $(\Sigma_{ij} + \Sigma_{ji})/2$. In a final step, the non-diagonal elements $\Sigma_{ij}$ were forced to values within the range $\pm\sqrt{\Sigma_{ii}\Sigma_{jj}}(1 - 10^{-5})$ to respect the properties of a covariance matrix. These steps did not cause any of the covariance matrices to become singular again.

The regularization of the covariance matrix is necessary to ensure that the marginalization integrals solved in Section 3.3

<footnote n="17">¹⁷ Here we use the term “regularized” in the sense where we ensure that the matrix is invertible, i.e., that it has a positive non-zero determinant.</footnote>
converge. Ill-defined covariance matrices with a negative determinant would cause the analytical solution to diverge.

The resulting moving group models are displayed in Figure 2, and their parameters are listed in Tables 1 and 9. The off-diagonal elements of the covariance matrices are provided in a FITS file with the BANYAN Σ algorithm (Gagné 2018a, 2018b). A kernel density estimate distribution was built for the 6D distribution of bona fide members using Silverman's rule of thumb, i.e., each data point is represented with a zero-covariance 6D multivariate Gaussian where the diagonal elements of the covariance matrix are given by:

$$ \Sigma_{ii} = (2N_k)^{-1/5}\sigma_i^2, \\ \Sigma_{ij} = 0, \quad i \neq j, \qquad (21) $$

where $\sigma_i$ is the standard deviation of the dimension i of the members' positions, and $N_k$ is the total number of members. The 1D projections of the kernel density estimate distribution are shown in the panels of Figure 2 that display histograms of the members' positions, and the residual difference between the multivariate Gaussian models and the 2D projections of the kernel density estimate distributions are shown in blue- and green-shaded backgrounds with the 2D distributions of members. There are several cases (such as TWA) where a number of projections show over- and under-densities of members by up to ≈40% of the multivariate Gaussian model, but we recommend against using multivariate Gaussian mixture models that would more correctly reproduce the distribution of known members until a large number of members are known (e.g., with the release of Gaia-DR2). Modeling the currently significantly incomplete distributions of association members with more complex models would negatively affect the ability of BANYAN Σ to recover the missing members that are located between the clumps of currently known members (i.e., in the blue-shaded regions in Figure 2).

The distance and radial velocity distributions of young moving group models were built by drawing 10⁵ synthetic objects from their spatial multivariate Gaussian model. The average distance or radial velocity was taken as the peak value of the resulting distribution and asymmetric characteristic widths covering half of 68% of the area under the curve on each side were measured. The resulting distance and radial velocity distributions as a function of the association ages are listed in Table 1 and displayed in Figures 3(a) and (b). These figures illustrate the range in distances and radial velocities where new moving group members of a given age can likely be discovered using BANYAN Σ.

In Figures 3(c) and (d), the characteristic spatial and kinematic scales $S_{\text{spa}}$ and $S_{\text{kin}}$ are displayed for different associations as a function of age, with:

$$ S_{\text{spa}} = \sqrt{|\Sigma_{\text{spa}}|^{1/3}}, \\ S_{\text{kin}} = \sqrt{|\Sigma_{\text{kin}}|^{1/3}}, $$

where $\Sigma_{\text{spa}}$ and $\Sigma_{\text{kin}}$ are 3 × 3 matrices that contain only the purely spatial or kinematic terms of the covariance matrix, respectively. The values of $S_{\text{spa}}$ and $S_{\text{kin}}$ are listed for each young association in Table 1.

Figure 3(c) illustrates how clusters are spatially much smaller than other types of associations, which become more dispersed as they age. The velocity dispersion of associations included in BANYAN Σ, as displayed in Figure 3(d), does not show a clear correlation with age. TAU is a clear outlier in both figures because a single multivariate Gaussian model is used to represent all of the TAU sub-groups.

The distribution of associations in the Galactic plane is displayed in Figures 4 and 5. The spatial location of the associations are defined as the contour that encompasses 68% of the projected multivariate Gaussian model, with an arbitrary minimum minor axis set at 4 pc for display. This figure illustrates how several associations are spatially close to each other, and how some of the nearest ones (3PMG, ABDMG) encompass the Sun. OCT is spatially the largest association because its members are spatially distributed in two distinct clumps, even though they share the same kinematics. Here we leave OCT as a single group because this may allow for the identification of new OCT members located spatially between the two clumps of known members.

## 6. A Model of Field Stars in the Solar Neighborhood

This section describes the construction of a kinematic model for the field hypothesis. It is based on the Besançon model (Robin et al. 1996, 2003, 2012, 2014, 2017)<fnref n="18" /> of the Galactic disk in the solar neighborhood (with a very small contribution from the Galactic halo), and uses the multivariate Gaussian formalism described in Section 3.1 for it to be compatible with the solution of the marginalization integrals developed in Section 3.3.

The Besançon Galactic model version used here follows the scheme described by Czekaj et al. (2014) for the thin disk population, which is based on their Model B (see their Table 5). In summary, thin disk stars are generated from a 3-slope initial mass function and a decreasing star formation rate, and follow the evolutionary tracks of Bertelli et al. (1994, 2008, 2009) for masses larger than $0.7M_{\odot}$, and from Chabrier & Baraffe (1997) for lower masses. Companions in binary systems are generated with a probability function that depends on the spectral type of the primary, and follow empirical mass-ratio and semimajor axis distributions, as described by Arenou (2011). The thick disk and halo populations are simulated with the best-fitting parameters obtained in the analysis of Robin et al. (2014) based on SDSS (Alam et al. 2015) and 2MASS data, and the isochrones of Bergbusch & Vandenberg (1992).

There are two complications that prevent a correct modeling of the field star kinematics with a simple multivariate Gaussian model: (1) the distribution in Z is similar to a hyperbolic secant function, which has wider wings than a Gaussian distribution; and (2) the distributions in X and Y are approximately uniform in the solar neighborhood.

The first problem can be addressed by modeling the field hypothesis with a mixture of N multivariate Gaussian

<footnote n="18">¹⁸ Available at http://model2016.obs-besancon.fr/.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Distance, radial velocity, spatial size and kinematic scatter of BANYAN Σ association models as a function of age. Associations of the solar neighborhood provide individual epochs that cover a large period of ages, relevant to disk evolution, planetary formation, and brown dwarfs' atmospheric cooling. Associations in our sample seem to display an increasing spatial size as a function of age, except for the denser open clusters. TAU is also an exception because its model includes several sub-groups. See Section 5 for more details.</figcaption></figure>

distributions:

$$ \mathcal{P}_{\text{field}} = \sum_{j=1}^{N} c_j \mathcal{P}_{j,\text{field}}, $$

$$ \mathcal{P}_{j,\text{field}} = \frac{1}{\sqrt{(2\pi)^6 |\bar{\Sigma}_j|}} \exp\left(-\frac{1}{2}(\bar{Q} - \tau_j)^T \bar{\Sigma}_j^{-1} (\bar{Q} - \bar{\tau}_j)\right), \\ \sum_{j=1}^{N} c_j = 1, $$

which yields the solution described in Equation (12) for the probability $\mathcal{P}({O_i}|H_{j,\text{field}})$ associated with field component j. The resulting field probability will then be:

$$ \mathcal{P}({O_i}|H_{\text{field}}) = \sum_{j=1}^{N} c_j \mathcal{P}({O_i}|H_{j,\text{field}}). $$

The second problem of the approximately uniform X and Y distributions can be mitigated by artificially inflating the two diagonal elements of all covariance matrices $\bar{\Sigma}_j$ corresponding to the X and Y dimensions by a factor much larger than the typical distances which will be involved in using BANYAN Σ. The density of the field model will however need to be re-adjusted to avoid affecting the stellar density in the solar neighborhood.

The very small covariance between all XYZUVW coordinates of field stars compared to their variances (the Pearson correlation coefficient of all dimensions is smaller than 0.1) makes the problem of fitting a mixture of multivariate Gaussians significantly easier. The covariance matrices $\bar{\Sigma}_j$ can be assumed diagonal and the fitting can be done simultaneously in four one-dimensional spaces Z, U, V and W instead of a single four-dimensional space. The X and Y
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Distribution in Galactic coordinates X and Y of all moving group models constructed in Section 3.1. The models included in BANYAN Σ cover all known associations and star-forming regions within 150 pc. This figure is an update of Figure 8 in Rice et al. (2011), although it is limited to 150 pc instead of 200 pc. See Section 5 for more details.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Distribution in Galactic coordinates X and Z of all moving group models constructed in Section 3.1. The color and linestyle coding is the same as that of Figure 4. See Section 5 for more details.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Unidimensional projections of the best-fitting ten-component multivariate Gaussian mixture models for the field stars (red dashed lines) compared to the distribution of stars in the Besançon Galactic model (black solid lines). The individual components of the best-fitting models are displayed as green lines, and best-fitting models that include fewer than seven Gaussian components are displayed as blue dashed lines. Space velocity distributions are slightly asymmetric and the distribution in the Z component of the Galactic position has much wider wings than a Gaussian model. The distributions in X and Y are approximated as uniform in the solar neighborhood. See Section 6 for more details.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Reduced $\chi^2$ of the multivariate field model as a function of the number of included multivariate Gaussian components. The distribution in space velocities UVW are slightly asymmetric and require $\gtrsim 7$ Gaussian components to be properly modeled, but the distribution in Z across the Galactic disk require more components as its wings are much wider than a Gaussian distribution. The distribution of stars in X and Y positions within the Galactic disk are approximated as uniform in the solar neighborhood. See Section 6 for more details.</figcaption></figure>

coordinates are ignored in a first step, as they will be approximated as uniform by using large Gaussian widths.

Multivariate Gaussians mixtures with $N = 1$ to 10 components were fitted to the Z, U, V, and W distributions of field stars using a Levenberg-Marquardt least-squares fit. The best-fitting models as well as the individual components of the $N = 10$ model are shown in Figure 6. Models with $N > 6$ provide a good visual fit to all Z, U, V, and W components.

In Figure 7, the reduced $\chi^2$ as a function of the number of mixture components $N$ is displayed for each dimension, and for the global fit across ZUVW. This figure shows that the goodness-of-fit does not improve significantly at $N > 7$ for the kinematic dimensions UVW, but Z keeps improving up to $N = 9–10$. The $N = 10$ components model was adopted, as it represents a good balance between accuracy and usability.

This simpler approach was preferred to one based on the Bayesian information criterion, as the very large number of field stars would have allowed for an arbitrarily large number of mixture components, which would make BANYAN $\Sigma$ impractical to use while causing very little difference in the calculated probabilities. The simplicity of the least-squares fitting also allowed the identification of a good general solution without needing to compute resources-intensive likelihood functions based on a very large number of field stars.
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Largest distance encompassing 90% of randomly selected field stars, as a function of lower cuts on total proper motion *μ* and absolute Galactic latitude |*b*|. The criteria used in the BASS survey for young brown dwarfs (Gagné et al. 2015b, 2015c) and the search for red L dwarfs of Schneider et al. (2017) are displayed with the star and triangle symbols, respectively. This figure demonstrates that field interlopers with distances up to ≈750 pc and ≈500 pc are likely contaminating the two respective input samples. BANYAN II has a limited capability of capturing contaminants at distances further than 200 pc, as more distant stars were not included in its model of the Galactic field. See Section 6 for more details.</figcaption></figure>

The *X* and *Y* components of the field model were arbitrarily set to a large characteristic width of 1500 pc. To ensure that this did not affect the density of stars in the solar neighborhood, a Monte Carlo simulation was used to draw field objects until a ratio of stars within 300 pc to the total number of stars could be counted with a precision of less than 1%, assuming Poisson error bars. This required a total of 10⁹ synthetic stars to be drawn and yielded a ratio of 1.69 × 10⁻⁵, which was divided to the total number of objects within 300 pc in the Besançon model (7.15 × 10⁶) to re-normalize the field model. This ensures that the multivariate Gaussian mixture model has the same density of stars as the Besançon model within 300 pc. The adopted field model parameters are provided as a FITS file containing all input data used by the BANYAN Σ IDL and Python codes (Gagné 2018a, 2018b).

## 6.1. The Spatial Size of Proper Motion and Galactic Latitude-limited Stellar Samples

One common way to eliminate distant stars from a sample is to impose a lower limit on the total proper motion and/or Galactic latitude. The model of the nearby Galactic disk developed in this work was used to determine the efficiency of such proper motion and Galactic latitude cuts at selecting nearby stars.

A set of 200 thresholds on the magnitude of proper motion and five thresholds on the Galactic latitude were selected uniformly in the ranges 5–300 mas yr⁻¹ and 0°–40°, respectively. For each combination of thresholds, the ZUVW coordinates of 10⁷ stars were drawn randomly from the multivariate Gaussians mixture model of the Galactic disk. Because the *X* and *Y* coordinates are approximated as uniform

in the solar neighborhood, they were drawn from a uniform random distribution bounded within a distance that produces a good sampling of the distance distribution of stars selected by the proper motion threshold. Bounds of ±10,000 pc/μ (mas yr⁻¹) on both *X* and *Y* were found to be adequate. All stars with proper motions and Galactic latitudes larger than the specific set of thresholds were selected, and the smallest distance that encompasses 90% of the sample was calculated.

The resulting distances encompassing 90% of a sample are displayed in Figure 8 as a function of proper motion and Galactic latitude selection cuts. This figure demonstrates that agressive cuts on proper motion (μ > 100 mas yr⁻¹) must be used to limit a sample to distances <~500 pc. The threshold on proper motion and Galactic latitude of the BANYAN All-Sky Survey for members of young associations (Gagné et al. 2015c) only limited their sample to distances of <800 pc, much larger than the 200 pc distance limit that was used to build the model of field stars in BANYAN II (Gagné et al. 2014). This problem was mitigated by the fact that the survey focused on substellar objects well detected in 2MASS, which limited the sample to distances <~200 pc for >M5-type objects. However, this outlines that a model of field stars that remains valid at much larger distances, such as the one developed in this section, is necessary to limit the rate of false-positives in all-sky searches for stellar members of young associations.

## 7. The Choice of Bayesian Priors

The contamination and recovery rates in a sample of candidate members selected with BANYAN Σ will be dependent on the young association where a given star is classified as a likely member. The more distant associations will provide a larger recovery rate of true association members at a fixed rate of contamination, because the members are distributed on a smaller region of the sky. As a consequence, using the same Bayesian probability threshold for the candidate members of all young associations will result in samples of vastly different sizes, completion, and contamination rates.

It is possible to adjust the Bayesian priors of the young associations in a way that equalizes the contamination or recovery rates of all young associations at an arbitrary Bayesian probability threshold. These priors will be used in the Bayesian membership probability determination (see Equation (1)). The threshold *P* = 90% was selected here so that BANYAN Σ approaches recovery rates *R*<sub>*k*</sub> of 50% (μ only), 68% (μ + ν), 82% (μ + ϖ), or 90% (μ + ν + ϖ) in terms of the fraction of recovered bona fide members. This decision is arbitrary, but will allow translating the BANYAN Σ probabilities to survey completeness fractions more simply. This will therefore make BANYAN Σ easier to use in searches for new members across all 27 associations.

This was done for each young association $H_k$ with a Monte Carlo method. The XYZUVW coordinates of 10⁷ synthetic stars were drawn from the kinematic model of the association. All coordinates were transformed to sky position, proper motion, radial velocity, and distance. Gaussian random error bars of 10 mas yr⁻¹ were added to each component of the proper motion, and radial velocity and distance were assumed to be missing. Simplified BANYAN Σ probabilities $P'(H_k|\{O_i\})$ were determined for all synthetic objects by ignoring all other
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Panel (a): receiver operating characteristic curves of different membership classification tools for distinguishing members of βPMG, THA, COL, TWA, and ABDMG from field objects based on only sky position and proper motion, and ignoring cross-contamination between associations. CP designates the convergent point tool. BANYAN Σ achieves the best performance under this metric, especially in the region TPR > 0.6. Individual values of the probability or goodness-of-fit metric thresholds are indicated along ROC curves. Panel (b): receiver operating characteristic curves of different membership classification tools for distinguishing members of the same five associations, using only sky position and proper motion, and ignoring field contamination. Classification tools that use more complex kinematic models of young associations such as BANYAN Σ tend to perform better under this metric. See Section 8 for more details.</figcaption></figure>

groups $H_j, j \neq k$:

$$P'(H_k|\{O_i\}) = \frac{P'(H_k)P'(\{O_i\}|H_k)}{P'(H_{\text{field}})P'(\{O_i\}|H_{\text{field}})},$$

where all initial priors $P'(H_k)$ and $P'(H_{\text{field}})$ were set to unity.

A set of 10⁵ probability thresholds $P_t$ were defined to explore how they affected the number of true positives ($N_{TP}$; association members with $P'(H_k|\{O_i\}) \ge P_t$) and false negatives ($N_{FN}$; association members with $P'(H_k|\{O_i\}) < P_t$).

The probability thresholds $P_t$ were distributed along:

$$P_t = \frac{1}{2} \left( \frac{\arctan(\omega x)}{\arctan(\omega)} + 1 \right), \\ \omega = 95,$$

where $x$ is an array of 10⁵ uniformly distributed values in the range [-1, 1]. This produces an array of thresholds where $P_t \approx 0$ and $P_t \approx 1$ are especially well sampled as $\omega$ takes larger values.

Recovery rates $R_k$ (also called “true positive rates” or TPRs) are determined as a function of threshold $P_t$, with:

$$R_k = \frac{N_{TP}}{N_{TP} + N_{FN}},$$

where $N_{TP}$ is the number of synthetic stars originating from the model of association association $H_k$ with $P_k \ge P_t$, and $N_{FN}$ is the number with $P_k < P_t$.

For each young association, the probability threshold $P_{t,crit}$ that generates the desired recovery rate (50%–90% depending on observables) was then selected, and a multiplicative factor $\alpha_k$ that ensures $P_{t,crit} = 90\%$ when included to the Bayesian

prior is then determined:

$$\alpha_k = \frac{0.9(1 - P_{t,crit})}{(1 - 0.9)P_{t,crit}}.$$

To avoid biasing the young association probabilities of ambiguous candidate members, an average of each young association factor $\alpha_k$, weighted by the individual membership probabilities $P(H_k|\{O_i\})$, is used to determine the field prior:

$$P(H_k) = 1, \\ P(H_{\text{field}}) = \left( \frac{\sum'_{k} \alpha_k P(H_k|\{O_i\})}{\sum'_{k} P(H_k|\{O_i\})} \right)^{-1},$$

where the sum $\sum'_k$ excludes $k = \text{field}$. The quantities $\alpha_k$ are fixed for all candidate stars, but $P(H_{\text{field}})$ has to be computed for each star, since it depends on $P(H_k|\{O_i\})$.

This choice of priors ensures that the recovery rates are similar (between 50% and 90% depending on the available observables) across all young associations when a probability threshold $P = 90\%$ is adopted, without biasing the relative young association membership probabilities. However, the false-positive rates at $P = 90\%$ will be different for each association, and are discussed in Section 8. The resulting $\alpha_k$ values are listed in Table 1.

## 8. The Performance of Banyan Σ as a Bayesian Classifier

This section describes the classification performance of BANYAN Σ (in normal mode and in UVW-only mode), and compares it to those of the convergent point tool (Jones 1971; de Bruijne 1999; Mamajek 2005; Rodriguez et al. 2011), BANYAN I (Malo et al. 2013), BANYAN II (Gagné et al. 2014), LACEwING (Riedel et al. 2017b) and the M-value metric introduced by Bowler et al. (2017). The probabilities or goodness-of-fit metrics of these different tools cannot be
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Left panel: false-positive rates for all young associations considered here as a function of input observables, for a threshold $P_t = 90\%$. The associations are sorted with decreasing false-positive rates in the case where only sky position and proper motion are considered. The grayed out region corresponds to the limit of our false-positive rate estimations, and corresponds to no false-positives in the full Besançon Galactic model described in Section 8. Right panel: false-positive rate as a function of the characteristic angular size of young associations, defined as the characteristic size of each young association $S_{spa}$ (see Section 5, Figure 3(c) and Table 1) divided by its distance. The false-positive rates displayed in this figure correspond to the case where only sky position and proper motion are considered in BANYAN $\Sigma$. There is a clear trend where associations with larger characteristic angular sizes suffer from more contamination from field stars. This effect alone dominates the false-positive rates, with only two exceptions: HYA and UMA are outliers with far fewer false-positives than would be expected from this trend. This is caused by their average UVW kinematics that differ significantly from those of most field stars and other young associations. See Section 8 for more details.</figcaption></figure>

compared in absolute terms; however their rate of contamination as a function of the rate of true member recovery are physically meaningful and can be directly compared.

The 7.15 × 10⁶ objects within 300 pc from the Besançon Galactic model and the bona fide members listed in Table 5 were used to determine the classification performance of each tool. In each case, the membership probability that each object belongs to a group $H_k$ or the field $H_{field}$ were calculated using only sky position and proper motion, and the number of recovered true members $N_{k,TP}$ from $H_k$ was counted for a range of thresholds $P_t$, with:

$$ \frac{P_k}{P_k + P_{field}} \geq P_t, \qquad (22) $$

or $P_k \geq P_t$ for the tools that do not have a field hypothesis.

The number of false positives $N_{k,FP}$ was defined as the number of Besançon objects that respect the same criterion. This particular way of normalizing probabilities ignores cross-contamination between young groups, and instead focuses on the contamination from field stars. These quantities make it possible to build receiver operating characteristic (ROC) curves, defined as the true positives rate TPR = $N_{k,TP}/N_k$ as a function of the false positives rate FPR = $N_{k,FP}/N_{field}$. The straight line defined by TPR = FPR corresponds to the performance of a random classification, and a ROC curve that is farthest above TPR > FPR corresponds to an optimal classification performance.

The ROC curve of each classification tool was built by taking the sum of $N_{k,TP}$ for only the young associations that are common to all tools, interpolated on a fixed array of $N_{k,FP}$. These young associations considered by all tools are βPMG, THA, COL, TWA, and ABDMG. The resulting ROC curves are displayed in Figure 9(a). The ROC curves do not compare

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Receiver operating characteristic curve for BANYAN $\Sigma$ as a function of the observables used, for all groups included here. Sky coordinates were used in all cases, and cross-contamination between the groups was ignored. The addition of radial velocity and distance cut down the rate of field contamination by factors of ~10 and ~100 respectively at a fixed recovery rate. See Section 8 for more details.</figcaption></figure>

the particular thresholds of each different tool, which are defined in different ways, but rather compares their astrophysically meaningful TPRs as a function of FPRs, which can be directly compared in an informative way. BANYAN $\Sigma$ achieves a performance slightly better than BANYAN II and BANYAN I, especially at large true-positive rates (TPR > 0.6). It is likely that the lack of a field model in LACEwING, the convergent point tool, and the M-value is the main reason they do not perform as well as the BANYAN tools under this metric. BANYAN $\Sigma$ in UVW-only mode
<figure><img src="image_14.png" /><figcaption>**Figure 12.** True positive rate (TPR) as a function of Bayesian probability threshold $P_t$ for different young associations, when only sky position and proper motion are considered. Bayesian priors were chosen for $P_t$ = 90% to yield TPR = 50% in this scenario (e.g., see the thick gray lines). The dashed thick gray lines represent the range in $P_t$ for which the range of recovery rates are reported in Table 10. See Section 8 for more details.</figcaption></figure>

achieves a much lower performance under this metric, but performs better than the *M*-value and the convergent point tool.

The resulting FPRs of all young associations are displayed in Figure 10(a) for each configuration of BANYAN Σ. In Figure 10(b), the FPRs in the case where only sky position and proper motion are used are displayed as a function of the characteristic angular size of the young associations, defined as their characteristic spatial size $S_{\text{spa}}$ (see Section 5) divided by their distance. The FPRs are dominated by the characteristic angular size of young associations: the nearby and large associations that cover a significant area of the sky are much harder to distinguish from field stars, because there is a much larger set of field stars that can match their kinematics by pure chance. HYA and UMA suffer from much less contamination that would be expected given their characteristic angular size; this is due to their average kinematics that significantly differ from most field stars and from other young associations.

Another set of ROC curves was built in a similar way to measure the cross-contamination performance, by ignoring the field hypothesis and defining false-positives as members recovered in association $H_k$ that originate from another association $H_l$, $l \neq k$. The resulting ROC curves are displayed in Figure 9(b). The classification tools that use more complex kinematic models tend to perform better under this metric: BANYAN Σ achieves the best performance, followed by BANYAN II, BANYAN I, the convergent point tool, BANYAN Σ in UVW-only mode, the *M*-value, and LACEwING. The kinematic models of LACEwING are similar to those of BANYAN II, and its lower performance may instead be related to approximations that are done when transforming the $N\sigma$ metrics taken in sky position and proper motion space to probabilities directly. In particular, most of the LACEwING cross-contamination is due to confusion between ABDMG and βPMG, the members of which span the widest distributions of sky positions and proper motions.

In order to characterize the classification gains that are obtained by adding more observables in BANYAN Σ, field contamination ROC curves were built for each mode: (1) proper motion only, (2) proper motion and radial velocity, (3) proper motion and distance, or (4) proper motion, radial velocity, and distance. These ROC curves were built from all associations available in BANYAN Σ, and are displayed in Figure 11. This figure demonstrates that using measurements of radial velocity and distance cut down field contamination by factors of ~10 and ~100 respectively at a fixed recovery rate. The fact that using only a distance measurement makes BANYAN Σ about a factor of ten better compared with only a radial velocity measurement is likely a consequence of distance being useful to constrain both XYZ and UVW sets of coordinates. Radial velocity only helps to constrain UVW. This observation was already made for BANYAN II and LACEwING (Gagné et al. 2014; Riedel et al. 2017b).

The models of young associations were also used to draw 10⁵ synthetic association members in order to obtain smooth TPRs as a function of probability threshold for each group in BANYAN Σ, which are displayed in Figure 12. This figure illustrates the effect of the $\alpha_k$ thresholds described in Section 7, causing the TPR curves of all young associations to meet at $P_t$ = 90% and TPR = 50%. Similar curves were built for FPRs and Matthews correlation coefficients (MCCs; Matthews & Czerwinski 1975), defined in the range [-1, 1] indicate the quality of a Bayesian classifier: An MCC of -1 indicates a perfect mis-classification, a value of 0 indicates a performance similar to a random classification, and a value of +1 indicates a perfect classifier. In these particular simulations, the number of true positives and false negatives (i.e., all synthetic objects originating from a young association model $H_k$) were scaled in the range [0, $N_k$] instead of [0, 10⁵], to give a realistic representations of the the FPR and MCC, where $N_k$ represents the number of stars in each young association. Because young associations with $N_k$ < 20 are expected to be incomplete, we have set a minimal value of $N_k$ = 20 in these situations. In a similar way, the number of true negatives and false positives (i.e., all synthetic objects originating from the Besançon Galactic model) were scaled in the range [0, 217,680], corresponding to the number of OBAFG-type stars in the model (out of 7.15 × 10⁶ stars). In doing this, we assume that most young associations are complete at this fraction; this approximation only affects the FPR and MCC values reported in this work. The values $N_k$ are listed for each young association in Table 1.

The resulting TPR, FPR, and MCC values for each young association, using each possible set of observables, are reported in Table 10, as the range of possible values within $P_t \in$ 85%–95% and centered on $P_t$ = 90%. They demonstrate how members of distant associations are much easier to distinguish from field objects because of their narrower distribution on the sky, with smaller FPRs and larger MCCs. Both the IDL and Python implementations of BANYAN Σ specify as an output for each star the TPR and FPR of a sample that would be constructed with only candidate members that have membership probabilities equal or higher than those of the star in question. This will allow users to interpret the Bayesian probabilities more easily without needing to rely on Table 10.

We investigated the effect of ignoring covariances between the input observables (see Section 3.3) with a 10⁴-element Monte Carlo simulation. The median proper motion and parallax errors of all bona fide members in Gaia-DR1
<table><caption>**Table 10**
BANYAN Σ Classification Performance for Different Young Associations as a Function of Input Observables</caption><thead><tr><th rowspan="2">Asso.</th><th colspan="4">TPR<sub>crit</sub></th><th colspan="4">log<sub>10</sub> FPR<sub>crit</sub></th><th colspan="4">log<sub>10</sub> MCC<sub>crit</sub></th></tr><tr><th>μ</th><th>μ, ν</th><th>μ, ρ</th><th>μ, ν, ρ</th><th>μ</th><th>μ, ν</th><th>μ, ρ</th><th>μ, ν, ρ</th><th>μ</th><th>μ, ν</th><th>μ, ρ</th><th>μ, ν, ρ</th></tr></thead><tbody><tr><td>118TAU</td><td>0.5<sup>+0.3</sup><sub>-0.1</sub></td><td>0.7<sup>+0.2</sup><sub>-0.1</sub></td><td>0.82<sup>+0.10</sup><sub>-0.05</sub></td><td>0.90<sup>+0.07</sup><sub>-0.03</sub></td><td>-4.9<sup>+0.3</sup><sub>-0.2</sub></td><td>-5.6 ± 0.2</td><td>-5.9<sup>+0.3</sup><sub>-0.2</sub></td><td>-6.4<sup>+0.5</sup><sub>-0.1</sub></td><td>-1.39<sup>+0.20</sup><sub>-0.01</sub></td><td>-0.88<sup>+0.08</sup><sub>-0.02</sub></td><td>-0.70<sup>+0.07</sup><sub>-0.06</sub></td><td>-0.43<sup>+0.10</sup><sub>-0.01</sub></td></tr><tr><td>ABDMG</td><td>0.5<sup>+0.2</sup><sub>-0.1</sub></td><td>0.68<sup>+0.10</sup><sub>-0.07</sub></td><td>0.82<sup>+0.20</sup><sub>-0.06</sub></td><td>0.90<sup>+0.09</sup><sub>-0.03</sub></td><td>-2.8<sup>+0.4</sup><sub>-0.2</sub></td><td>-3.5<sup>+0.3</sup><sub>-0.2</sub></td><td>-4.0<sup>+0.3</sup><sub>-0.2</sub></td><td>-4.5<sup>+0.3</sup><sub>-0.1</sub></td><td>-2.24<sup>-0.01</sup></td><td>-1.78<sup>+0.07</sup><sub>-0.05</sub></td><td>-1.43<sup>+0.07</sup><sub>-0.05</sub></td><td>-1.14<sup>+0.08</sup><sub>-0.06</sub></td></tr><tr><td>BPMG</td><td>0.5<sup>+0.2</sup><sub>-0.1</sub></td><td>0.68<sup>+0.10</sup><sub>-0.07</sub></td><td>0.82<sup>+0.10</sup><sub>-0.06</sub></td><td>0.90<sup>+0.08</sup><sub>-0.03</sub></td><td>-2.6<sup>+0.4</sup><sub>-0.2</sub></td><td>-3.2<sup>+0.3</sup><sub>-0.2</sub></td><td>-4.4<sup>+0.3</sup><sub>-0.1</sub></td><td>-4.8<sup>+0.3</sup><sub>-0.2</sub></td><td>-2.39<sup>-0.02</sup></td><td>-1.95<sup>+0.06</sup><sub>-0.04</sub></td><td>-1.28<sup>+0.06</sup><sub>-0.05</sub></td><td>-1.03 ± 0.09</td></tr><tr><td>CAR</td><td>0.5<sup>+0.2</sup><sub>-0.1</sub></td><td>0.68<sup>+0.20</sup><sub>-0.09</sub></td><td>0.82<sup>+0.10</sup><sub>-0.05</sub></td><td>0.90<sup>+0.07</sup><sub>-0.03</sub></td><td>-2.9<sup>+0.5</sup><sub>-0.2</sub></td><td>-3.6<sup>+0.4</sup><sub>-0.2</sub></td><td>-5.0<sup>+0.4</sup><sub>-0.2</sub></td><td>-5.3<sup>+0.5</sup><sub>-0.1</sub></td><td>-2.36<sup>+0.06</sup></td><td>-1.93<sup>+0.03</sup><sub>-0.04</sub></td><td>-1.12<sup>+0.10</sup><sub>-0.08</sub></td><td>-0.96<sup>+0.20</sup><sub>-0.04</sub></td></tr><tr><td>CARN</td><td>0.5<sup>+0.2</sup><sub>-0.1</sub></td><td>0.68<sup>+0.20</sup><sub>-0.08</sub></td><td>0.82<sup>+0.10</sup><sub>-0.06</sub></td><td>0.90<sup>+0.09</sup><sub>-0.03</sub></td><td>-3.3<sup>+0.5</sup><sub>-0.2</sub></td><td>-4.0<sup>+0.3</sup><sub>-0.2</sub></td><td>-4.5<sup>+0.3</sup><sub>-0.2</sub></td><td>-4.9<sup>+0.3</sup><sub>-0.1</sub></td><td>-2.18<sup>+0.03</sup></td><td>-1.71<sup>+0.03</sup><sub>-0.05</sub></td><td>-1.39<sup>+0.05</sup><sub>-0.06</sub></td><td>-1.14<sup>+0.10</sup><sub>-0.06</sub></td></tr><tr><td>CBER</td><td>0.5<sup>+0.5</sup><sub>-0.2</sub></td><td>0.7<sup>+0.4</sup><sub>-0.1</sub></td><td>0.82<sup>+0.20</sup><sub>-0.06</sub></td><td>0.90<sup>+0.10</sup><sub>-0.04</sub></td><td>-4.5<sup>+0.6</sup><sub>-0.2</sub></td><td>-5.2<sup>+0.6</sup><sub>-0.2</sub></td><td>-6.9<sup>-0.3</sup></td><td>&lt;-7</td><td>-1.40<sup>+0.80</sup><sub>-0.02</sub></td><td>-0.96<sup>+0.05</sup><sub>-0.03</sub></td><td>-0.19<sup>+0.08</sup><sub>-0.06</sub></td><td>-0.17<sup>+0.02</sup><sub>-0.01</sub></td></tr><tr><td>COL</td><td>0.50<sup>+0.10</sup><sub>-0.08</sub></td><td>0.68<sup>+0.09</sup><sub>-0.05</sub></td><td>0.82<sup>+0.10</sup><sub>-0.05</sub></td><td>0.90<sup>+0.06</sup><sub>-0.03</sub></td><td>-2.7<sup>+0.3</sup><sub>-0.1</sub></td><td>-3.2<sup>+0.2</sup><sub>-0.1</sub></td><td>-4.5<sup>+0.3</sup><sub>-0.1</sub></td><td>-4.8<sup>+0.2</sup><sub>-0.1</sub></td><td>-2.46<sup>-0.01</sup></td><td>-2.08<sup>+0.05</sup><sub>-0.03</sub></td><td>-1.33<sup>+0.07</sup><sub>-0.05</sub></td><td>-1.13<sup>+0.08</sup><sub>-0.05</sub></td></tr><tr><td>CRA</td><td>0.5<sup>+0.2</sup><sub>-0.1</sub></td><td>0.68<sup>+0.20</sup><sub>-0.09</sub></td><td>0.82<sup>+0.20</sup><sub>-0.06</sub></td><td>0.90<sup>+0.09</sup><sub>-0.03</sub></td><td>-5.06<sup>+0.30</sup><sub>-0.09</sub></td><td>-5.8 ± 0.1</td><td>-6.55<sup>-0.20</sup></td><td>&lt;-7</td><td>-1.31<sup>+0.10</sup><sub>-0.06</sub></td><td>-0.81<sup>+0.08</sup></td><td>-0.39<sup>+0.04</sup><sub>-0.05</sub></td><td>&gt;-0.22</td></tr><tr><td>EPSC</td><td>0.5<sup>+0.3</sup><sub>-0.2</sub></td><td>0.7<sup>+0.3</sup><sub>-0.1</sub></td><td>0.82<sup>+0.20</sup><sub>-0.06</sub></td><td>0.90<sup>+0.10</sup><sub>-0.03</sub></td><td>-4.9<sup>+0.6</sup><sub>-0.2</sub></td><td>-5.3<sup>+0.5</sup><sub>-0.2</sub></td><td>&lt;-7</td><td>-6.85<sup>-0.50</sup></td><td>-1.35<sup>+0.20</sup><sub>-0.01</sub></td><td>-1.00 ± 0.05</td><td>-0.43<sup>+0.10</sup><sub>-0.03</sub></td><td>-0.22<sup>+0.03</sup><sub>-0.20</sub></td></tr><tr><td>ETAC</td><td>0.5<sup>+0.3</sup><sub>-0.2</sub></td><td>0.7<sup>+0.3</sup><sub>-0.1</sub></td><td>0.82<sup>+0.20</sup><sub>-0.07</sub></td><td>0.90<sup>+0.10</sup><sub>-0.04</sub></td><td>-6.0<sup>+0.4</sup><sub>-0.2</sub></td><td>-6.37<sup>+0.01</sup><sub>-0.40</sub></td><td>&lt;-7</td><td>&lt;-7</td><td>-0.85<sup>+0.30</sup><sub>-0.03</sub></td><td>-0.6<sup>+0.2</sup><sub>-0.1</sub></td><td>&gt;-0.22</td><td>&gt;-0.22</td></tr><tr><td>HYA</td><td>0.5<sup>+0.3</sup><sub>-0.1</sub></td><td>0.68<sup>+0.20</sup><sub>-0.09</sub></td><td>0.82<sup>+0.10</sup><sub>-0.06</sub></td><td>0.90<sup>+0.09</sup><sub>-0.03</sub></td><td>-5.4<sup>+0.6</sup><sub>-0.2</sub></td><td>-6.2<sup>+0.7</sup><sub>-0.3</sub></td><td>-6.6<sup>+0.3</sup><sub>-0.2</sub></td><td>-6.82<sup>-0.40</sup></td><td>-0.67<sup>+0.05</sup><sub>-0.02</sub></td><td>-0.27<sup>+0.04</sup><sub>-0.06</sub></td><td>-0.12<sup>+0.02</sup><sub>-0.01</sub></td><td>-0.07<sup>+0.02</sup><sub>-0.05</sub></td></tr><tr><td>IC2391</td><td>0.5<sup>+0.3</sup><sub>-0.1</sub></td><td>0.68<sup>+0.20</sup><sub>-0.07</sub></td><td>0.82<sup>+0.20</sup><sub>-0.04</sub></td><td>0.90<sup>+0.05</sup><sub>-0.02</sub></td><td>-5.3<sup>+0.3</sup><sub>-0.1</sub></td><td>-5.7<sup>+0.3</sup><sub>-0.1</sub></td><td>-6.0<sup>+0.4</sup></td><td>-6.6<sup>+0.3</sup></td><td>-1.18<sup>+0.20</sup><sub>-0.03</sub></td><td>-0.88<sup>+0.03</sup></td><td>-0.66<sup>+0.10</sup><sub>-0.02</sub></td><td>-0.36<sup>+0.09</sup></td></tr><tr><td>IC2602</td><td>0.5<sup>+0.3</sup><sub>-0.1</sub></td><td>0.68<sup>+0.20</sup><sub>-0.08</sub></td><td>0.82<sup>+0.10</sup><sub>-0.05</sub></td><td>0.90<sup>+0.06</sup><sub>-0.02</sub></td><td>-5.0<sup>+0.3</sup><sub>-0.1</sub></td><td>-5.42<sup>+0.20</sup><sub>-0.07</sub></td><td>&lt;-7</td><td>&lt;-7</td><td>-1.35<sup>+0.30</sup><sub>-0.05</sub></td><td>-1.00<sup>+0.08</sup><sub>-0.01</sub></td><td>-0.28 ± 0.02</td><td>&gt;-0.24</td></tr><tr><td>LCC</td><td>0.5<sup>+0.3</sup><sub>-0.1</sub></td><td>0.68<sup>+0.20</sup><sub>-0.09</sub></td><td>0.82<sup>+0.10</sup><sub>-0.06</sub></td><td>0.90<sup>+0.08</sup><sub>-0.03</sub></td><td>-3.2<sup>+0.4</sup><sub>-0.2</sub></td><td>-3.5<sup>+0.3</sup><sub>-0.2</sub></td><td>-4.5<sup>+0.3</sup><sub>-0.2</sub></td><td>-4.6 ± 0.2</td><td>-1.91<sup>+0.10</sup><sub>-0.02</sub></td><td>-1.62<sup>+0.02</sup><sub>-0.03</sub></td><td>-1.06<sup>+0.07</sup><sub>-0.06</sub></td><td>-0.97 ± 0.07</td></tr><tr><td>OCT</td><td>0.50<sup>+0.10</sup><sub>-0.06</sub></td><td>0.68<sup>+0.05</sup><sub>-0.03</sub></td><td>0.82<sup>+0.04</sup><sub>-0.02</sub></td><td>0.90<sup>+0.02</sup><sub>-0.01</sub></td><td>-3.0<sup>+0.1</sup><sub>-0.1</sub></td><td>-3.22<sup>+0.20</sup><sub>-0.09</sub></td><td>-3.43<sup>+0.20</sup><sub>-0.09</sub></td><td>-3.38<sup>+0.10</sup><sub>-0.05</sub></td><td>-2.31<sup>+0.01</sup></td><td>-2.09<sup>+0.04</sup><sub>-0.03</sub></td><td>-1.91<sup>+0.06</sup><sub>-0.03</sub></td><td>-1.89<sup>+0.04</sup><sub>-0.02</sub></td></tr><tr><td>PL8</td><td>0.5<sup>+0.3</sup><sub>-0.1</sub></td><td>0.68<sup>+0.20</sup><sub>-0.08</sub></td><td>0.82<sup>+0.20</sup><sub>-0.04</sub></td><td>0.90<sup>+0.05</sup><sub>-0.02</sub></td><td>-4.3<sup>+0.2</sup><sub>-0.1</sub></td><td>-4.6<sup>+0.2</sup><sub>-0.1</sub></td><td>-5.4<sup>+0.2</sup><sub>-0.1</sub></td><td>-5.7<sup>+0.2</sup><sub>-0.1</sub></td><td>-1.70<sup>+0.30</sup><sub>-0.03</sub></td><td>-1.38<sup>+0.03</sup></td><td>-0.94 ± 0.05</td><td>-0.72<sup>+0.08</sup><sub>-0.05</sub></td></tr><tr><td>PLE</td><td>0.5<sup>+0.4</sup><sub>-0.2</sub></td><td>0.7<sup>+0.3</sup><sub>-0.1</sub></td><td>0.82<sup>+0.10</sup><sub>-0.05</sub></td><td>0.90<sup>+0.09</sup><sub>-0.03</sub></td><td>-5.3<sup>+0.4</sup><sub>-0.3</sub></td><td>-5.7<sup>+0.7</sup><sub>-0.2</sub></td><td>-6.0<sup>+0.2</sup><sub>-0.4</sub></td><td>-6.4<sup>+0.5</sup><sub>-0.3</sub></td><td>-0.72<sup>+0.40</sup><sub>-0.01</sub></td><td>-0.43<sup>+0.03</sup><sub>-0.01</sub></td><td>-0.2 ± 0.1</td><td>-0.12<sup>+0.03</sup><sub>-0.06</sub></td></tr><tr><td>ROPH</td><td>0.5<sup>+0.2</sup><sub>-0.1</sub></td><td>0.68<sup>+0.20</sup><sub>-0.07</sub></td><td>0.82<sup>+0.20</sup><sub>-0.07</sub></td><td>0.90<sup>+0.10</sup><sub>-0.04</sub></td><td>-5.4 ± 0.2</td><td>-6.3 ± 0.1</td><td>&lt;-7</td><td>&lt;-7</td><td>-0.68<sup>+0.08</sup></td><td>-0.24<sup>+0.05</sup></td><td>&gt;-0.05</td><td>&gt;-0.04</td></tr><tr><td>TAU</td><td>0.50<sup>+0.10</sup><sub>-0.08</sub></td><td>0.68<sup>+0.10</sup><sub>-0.07</sub></td><td>0.82<sup>+0.20</sup><sub>-0.06</sub></td><td>0.90<sup>+0.09</sup><sub>-0.03</sub></td><td>-2.5<sup>+0.2</sup><sub>-0.1</sub></td><td>-3.0<sup>+0.2</sup><sub>-0.1</sub></td><td>-4.2<sup>+0.3</sup><sub>-0.2</sub></td><td>-4.5<sup>+0.3</sup><sub>-0.2</sub></td><td>-2.18<sup>+0.05</sup><sub>-0.01</sub></td><td>&gt;-1.81</td><td>-1.14<sup>+0.07</sup><sub>-0.06</sub></td><td>-0.93<sup>+0.10</sup><sub>-0.06</sub></td></tr><tr><td>THA</td><td>0.5<sup>+0.2</sup><sub>-0.1</sub></td><td>0.68<sup>+0.20</sup><sub>-0.08</sub></td><td>0.82<sup>+0.10</sup><sub>-0.05</sub></td><td>0.90<sup>+0.07</sup><sub>-0.03</sub></td><td>-3.9<sup>+0.5</sup><sub>-0.2</sub></td><td>-4.5<sup>+0.4</sup><sub>-0.2</sub></td><td>-5.3<sup>+0.3</sup><sub>-0.1</sub></td><td>-5.47<sup>+0.20</sup><sub>-0.08</sub></td><td>-1.76<sup>-0.03</sup></td><td>-1.33<sup>+0.08</sup><sub>-0.06</sub></td><td>-0.84<sup>+0.07</sup><sub>-0.04</sub></td><td>-0.71<sup>+0.05</sup><sub>-0.03</sub></td></tr><tr><td>THOR</td><td>0.5<sup>+0.3</sup><sub>-0.1</sub></td><td>0.7<sup>+0.2</sup><sub>-0.1</sub></td><td>0.82<sup>+0.20</sup><sub>-0.06</sub></td><td>0.90<sup>+0.10</sup><sub>-0.03</sub></td><td>-3.6<sup>+0.5</sup><sub>-0.2</sub></td><td>-4.2<sup>+0.4</sup><sub>-0.2</sub></td><td>-5.7<sup>+0.5</sup><sub>-0.3</sub></td><td>-6.1<sup>+0.5</sup><sub>-0.2</sub></td><td>-1.90<sup>+0.10</sup><sub>-0.01</sub></td><td>-1.49<sup>+0.01</sup><sub>-0.03</sub></td><td>-0.66<sup>+0.10</sup><sub>-0.09</sub></td><td>-0.46<sup>+0.10</sup><sub>-0.08</sub></td></tr><tr><td>TWA</td><td>0.5<sup>+0.2</sup><sub>-0.1</sub></td><td>0.68<sup>+0.20</sup><sub>-0.09</sub></td><td>0.82<sup>+0.10</sup><sub>-0.05</sub></td><td>0.90<sup>+0.08</sup><sub>-0.03</sub></td><td>-4.1<sup>+0.5</sup><sub>-0.2</sub></td><td>-4.6<sup>+0.4</sup><sub>-0.2</sub></td><td>-5.5<sup>+0.3</sup><sub>-0.1</sub></td><td>-5.8<sup>+0.3</sup><sub>-0.2</sub></td><td>-1.76<sup>+0.03</sup></td><td>-1.37 ± 0.04</td><td>-0.86<sup>+0.07</sup><sub>-0.04</sub></td><td>-0.67<sup>+0.10</sup><sub>-0.08</sub></td></tr><tr><td>UCL</td><td>0.5<sup>+0.3</sup><sub>-0.1</sub></td><td>0.68<sup>+0.20</sup><sub>-0.09</sub></td><td>0.82<sup>+0.10</sup><sub>-0.06</sub></td><td>0.90<sup>+0.08</sup><sub>-0.03</sub></td><td>-2.9<sup>+0.4</sup><sub>-0.2</sub></td><td>-3.2<sup>+0.3</sup><sub>-0.2</sub></td><td>-3.7<sup>+0.3</sup><sub>-0.2</sub></td><td>-3.9<sup>+0.2</sup><sub>-0.1</sub></td><td>-2.05<sup>+0.10</sup><sub>-0.01</sub></td><td>-1.75<sup>-0.02</sup></td><td>-1.40<sup>+0.09</sup><sub>-0.05</sub></td><td>-1.27<sup>+0.08</sup><sub>-0.06</sub></td></tr><tr><td>UCRA</td><td>0.5<sup>+0.2</sup><sub>-0.1</sub></td><td>0.68<sup>+0.10</sup><sub>-0.06</sub></td><td>0.82<sup>+0.07</sup><sub>-0.03</sub></td><td>0.90<sup>+0.03</sup><sub>-0.02</sub></td><td>-4.3<sup>+0.3</sup><sub>-0.1</sub></td><td>-4.8<sup>+0.2</sup><sub>-0.1</sub></td><td>-5.2<sup>+0.2</sup><sub>-0.1</sub></td><td>-5.32<sup>+0.20</sup><sub>-0.07</sub></td><td>-1.67<sup>+0.10</sup><sub>-0.02</sub></td><td>-1.29<sup>-0.02</sup></td><td>-1.02<sup>+0.07</sup><sub>-0.05</sub></td><td>-0.92<sup>+0.09</sup><sub>-0.03</sub></td></tr><tr><td>UMA</td><td>0.5<sup>+0.2</sup><sub>-0.1</sub></td><td>0.68<sup>+0.20</sup><sub>-0.09</sub></td><td>0.82<sup>+0.10</sup><sub>-0.06</sub></td><td>0.90<sup>+0.09</sup><sub>-0.03</sub></td><td>&lt;-7</td><td>&lt;-7</td><td>&lt;-7</td><td>&lt;-7</td><td>-0.58<sup>+0.09</sup><sub>-0.07</sub></td><td>-0.33 ± 0.02</td><td>&gt;-0.21</td><td>~0</td></tr><tr><td>USCO</td><td>0.5<sup>+0.2</sup><sub>-0.1</sub></td><td>0.68<sup>+0.20</sup><sub>-0.09</sub></td><td>0.82<sup>+0.10</sup><sub>-0.05</sub></td><td>0.90<sup>+0.08</sup><sub>-0.03</sub></td><td>-3.7<sup>+0.3</sup><sub>-0.2</sub></td><td>-4.0<sup>+0.3</sup><sub>-0.1</sub></td><td>-4.5<sup>+0.3</sup><sub>-0.1</sub></td><td>-4.7<sup>+0.3</sup><sub>-0.1</sub></td><td>-1.70<sup>+0.10</sup><sub>-0.02</sub></td><td>-1.39<sup>+0.02</sup><sub>-0.01</sub></td><td>-1.04<sup>+0.08</sup><sub>-0.04</sub></td><td>-0.93<sup>+0.09</sup><sub>-0.04</sub></td></tr><tr><td>XFOR</td><td>0.5<sup>+0.3</sup><sub>-0.1</sub></td><td>0.68<sup>+0.20</sup><sub>-0.08</sub></td><td>0.82<sup>+0.10</sup><sub>-0.04</sub></td><td>0.90<sup>+0.05</sup><sub>-0.02</sub></td><td>-5.4<sup>+0.4</sup><sub>-0.1</sub></td><td>-5.9<sup>+0.2</sup><sub>-0.3</sub></td><td>-6.82<sup>-0.30</sup></td><td>&lt;-7</td><td>-1.15<sup>+0.20</sup><sub>-0.03</sub></td><td>-0.78<sup>+0.06</sup><sub>-0.08</sub></td><td>-0.29<sup>+0.02</sup><sub>-0.08</sub></td><td>-0.25<sup>+0.02</sup></td></tr></tbody></table>

Note. True-positive rates (TPRs), false-positive rates (FPRs) and Matthews correlation coefficients (MCCs) are reported at the “critical” $P_t = 90\%$ threshold, and their reported ranges are for thresholds in the range 85%–95%. See Section 8 for more details.
<table><caption>**Table 11**
BANYAN Σ Probabilities for New Candidate Members of Sco-Cen Identified by Rizzuto et al. (2011)</caption><thead><tr><th>2MASS<br>Designation</th><th>Asso.</th><th>P<br>(%)</th></tr></thead><tbody><tr><td colspan="3"><em>Unambiguous candidate members of USCO</em></td></tr><tr><td>15480330-2512562</td><td>USCO</td><td>96.8</td></tr><tr><td>15574880-2331383</td><td>USCO</td><td>99.0</td></tr><tr><td>16052655-1948066</td><td>USCO</td><td>99.2</td></tr><tr><td>16120593-2314445</td><td>USCO</td><td>98.6</td></tr><tr><td>16203056-2006518</td><td>USCO</td><td>99.2</td></tr><tr><td>16301246-2506548</td><td>USCO</td><td>97.1</td></tr><tr><td colspan="3"><em>Unambiguous candidate members of UCL</em></td></tr><tr><td>13514960-3259387</td><td>UCL</td><td>96.3</td></tr><tr><td>14170338-3432122</td><td>UCL</td><td>96.6</td></tr><tr><td>14250102-3726493</td><td>UCL</td><td>97.8</td></tr><tr><td>14281043-2929299</td><td>UCL</td><td>97.3</td></tr><tr><td>14353043-4209281</td><td>UCL</td><td>99.6</td></tr><tr><td>14353149-4131026</td><td>UCL</td><td>99.9</td></tr><tr><td>15025928-3238357</td><td>UCL</td><td>96.9</td></tr><tr><td>15130106-3714480</td><td>UCL</td><td>99.6</td></tr><tr><td>15183199-4752307</td><td>UCL</td><td>97.6</td></tr><tr><td>15195653-3006249</td><td>UCL</td><td>95.4</td></tr><tr><td>15342085-3920572</td><td>UCL</td><td>99.6</td></tr><tr><td>15383263-3909384</td><td>UCL</td><td>99.8</td></tr><tr><td>15500707-5312351</td><td>UCL</td><td>96.0</td></tr><tr><td>16003131-3605164</td><td>UCL</td><td>98.6</td></tr><tr><td>16040712-3510367</td><td>UCL</td><td>96.8</td></tr><tr><td>16312294-3442153</td><td>UCL</td><td>95.6</td></tr><tr><td>16383094-3909083</td><td>UCL</td><td>95.5</td></tr><tr><td colspan="3"><em>Unambiguous candidate members of LCC</em></td></tr><tr><td>11454479-5241258</td><td>LCC</td><td>95.1</td></tr><tr><td>12044525-5915117</td><td>LCC</td><td>99.8</td></tr><tr><td>12113912-5222065</td><td>LCC</td><td>98.4</td></tr><tr><td>12263615-6305571</td><td>LCC</td><td>99.8</td></tr><tr><td>12474326-5941194</td><td>LCC</td><td>99.7</td></tr><tr><td colspan="3"><em>Ambiguous candidate members</em></td></tr><tr><td>13172895-4255587</td><td>TWA,UCL,LCC</td><td>79.4,10.5,5.7</td></tr><tr><td>13324248-5549391</td><td>LCC,UCL</td><td>89.4,7.8</td></tr><tr><td>13414477-5433339</td><td>LCC,UCL</td><td>76.6,21.6</td></tr><tr><td>14074081-4842144</td><td>UCL,LCC</td><td>80.4,18.8</td></tr><tr><td>15585013-3203082</td><td>UCL,USCO</td><td>74.3,23.4</td></tr><tr><td>16115069-2733098</td><td>USCO,UCL</td><td>49.3,47.4</td></tr><tr><td>16201552-2843008</td><td>USCO,UCL</td><td>77.7,20.6</td></tr></tbody></table>

Note. See Section 9 for more details.

(respectively 0.2 and 0.15 mas yr⁻¹, and 0.3 mas) were assigned as the measurement errors to the observables of a typical member (TWA 1). The radial velocity of TWA 1 was ignored to maximize the effect of the covariances in the other measurements, and therefore to be as conservative as possible. Random measurements of its proper motion and parallax were taken from a 3D multivariate Gaussian distribution, where each dimension corresponds to the two components of proper motion and the parallax, and assuming a 99.9% correlation between all dimensions. Each of these synthetic stars was assigned a vanishingly small error on its proper motion and distance (which we set to the median Gaia-DR1 values divided by one hundred). The membership probabilities of the 10⁴ synthetic stars were then calculated with BANYAN Σ, and the average probability was calculated, which corresponds to the final probability marginalized over all values of proper motion

and parallax. Ignoring the covariances in the observed measurements resulted in a negligible difference of ~0.09% in the membership probability, and similarly negligible differences of respectively 0.2% and 0.01% in the measured optimal distance and radial velocity. A similar Monte Carlo analysis was performed where the distance measurement was ignored, and yielded even smaller differences: 3 × 10⁻⁵% in probability, 3 × 10⁻⁶% in optimal distance and 8 × 10⁻⁶% in optimal radial velocity. The effects of ignoring covariances between the components of sky position and other quantities would be even smaller, because the error bars on the sky position are always significantly smaller than those on proper motion and parallax.

## 9. Classifying Previously Ambiguous Members with Banyan Σ

In this section, BANYAN Σ is used to assign classifications to previously ambiguous members encountered in the construction of the bona fide members lists (e.g., see the discussion in Section 4). None of the stars analyzed in this section was used to construct the young association models of BANYAN Σ. Only objects with a *P* ≥ 95% Bayesian probability of belonging to young associations are discussed here.

TWA 19 AB is a young star that was defined as a candidate member of TWA by Mamajek (2005), but later listed as a likely contaminant from LCC by Gagné et al. (2017b) based on a preliminary version of BANYAN Σ. Here we obtain a membership probability of *P* = 99% that it is a member of LCC.

Zuckerman & Song (2004) listed several tentative members of ABDMG as having a “questionable membership.” We find that five of them (HD 6569, HD 13482, HD 139751, HD 218860, and HD 224228) obtain a *P* > 97% membership probability associated with ABDMG.

2MASS J05361998–1920396 is a young L2 γ substellar object that was listed as an ambiguous member of COL and β PMG by Faherty et al. (2016). Here, we obtain an unambiguous *P* = 99.9% probability that it is a member of COL.

CP–68 1388, GSC 09235–01702, CD–69 1055, and MP Mus were all identified as ambiguous members of LCC or EPSC based on literature compilations. All but GSC 09235–01702 obtain an unambiguous *P* > 98% LCC membership probability when analyzed with BANYAN Σ. GSC 09235–01702 remains somewhat ambiguous with a 85% EPSC membership probability and a 12.7% LCC membership probability.

As mentioned in Section 4, there is a subset of new candidate members of the Sco-Cen region (consisting of the USCO, LCC and UCL sub-groups) discovered by Rizzuto et al. (2011) that were not assigned to either of its sub-groups. Of these candidate members, 35 benefit from at least a radial velocity and/or parallax measurement, and obtain a *P* > 95% young association membership probability. These objects are listed in Table 11 along with their respective probabilities in UCL, LCC and USCO.

Only seven of them remain ambiguous between more than one of the Sco-Cen subgroups, and they all obtain negligible membership probabilities for associations outside of the Sco-Cen region, with one exception: 2MASS J13172895–4255587 (F3 V; Houk 1978) obtains respective membership probabilities of 79.4% for TWA, 10.5% for UCL, 5.7% for LCC, and 4.3% for the field. This star has not been identified by Gagné et al. (2017b), who used a preliminary BANYAN Σ model of TWA
to identify new members based on *Hipparcos*. Further studies of this star may be able to confirm whether its age and radial velocity match those of TWA better than the Sco-Cen region.

## 10. Summary and Conclusions

A new Bayesian algorithm to identify young association members was presented. It derives membership probabilities from the sky position and proper motion of an object, and optionally radial velocity, parallax, and spectrophotometric distance constraints. It includes various improvements over its predecessor BANYAN II (Gagné et al. 2014), as it includes more associations, an updated list of bona fide members including Gaia-DR1 data, better spatial-kinematic models, a more accurate model of field stars, fewer approximations, new options, and a significantly enhanced execution speed due to an analytical solution of the marginalization integrals. One limitation of BANYAN Σ is that it cannot account for correlations in measured error bars, such as those reported in Gaia, but this results in biases of less than 0.1% in membership probability and measurements of the optimal distance and radial velocity.

The new BANYAN Σ tool includes all of the 27 currently known young associations within 150 pc, for which the current census of bona fide members is updated. An online tool is also made publicly available at http://www.exoplanetes.umontreal.ca/banyan/banyansigma.php. Additional figures and information on this work can be found on the website http://www.astro.umontreal.ca/~gagne.

We strongly encourage users to investigate the position of candidate members in appropriate color-magnitude diagrams using the BANYAN Σ optimal distance when no parallax measurements are used as inputs, as candidate lists generated without them suffer from ~100 times more false positives in comparison (see Section 8 and Figure 11). The BANYAN Σ algorithm can include distance constraints from population sequences in color-magnitude diagrams, and future versions will provide examples of such sequences.

This first version of BANYAN Σ will be the basis of a search for new isolated planetary-mass members in young associations based on 2MASS and AllWISE through the BASS-Ultracool survey (J. Gagné et al. 2018, in preparation; see also Gagné et al. 2015a, 2017a for preliminary results from BASS-Ultracool), as well as a search for new stellar members based on Gaia-DR1 (J. Gagné et al. 2018, in preparation). The release of Gaia-DR2 will allow us to significantly improve the spatial and kinematic models of BANYAN Σ, and identify new members of young associations that are not part of the Tycho catalog. A second version of the BANYAN Σ software with such improved models will be released in a future work that will aim at furthering the census of young association members based on Gaia-DR2.

The authors would like to thank the anonymous referee and the AAS statistics consultant for valuable and detailed comments that significantly improved the quality of this paper. We thank Bruno S. Alessi and Eric Bubar for sharing data. We thank Noé Aubin-Cadot, Joel Kastner, Thierry Bazier-Matte, Simon Gélinas, Jean-François Désilets and Brendan Bowler for useful comments, as well as the anonymous referee of Gagné et al. (2015b), who suggested the inclusion of a parallax motion correction in the BANYAN tools.

This research 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; Skrutskie et al. 2006), 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); data products from the Wide-field Infrared Survey Explorer (WISE; and Wright et al. 2010), which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory (JPL)/Caltech, funded by NASA. This project was developed in part at the 2017 Heidelberg Gaia Sprint, hosted by the Max-Planck-Institut für Astronomie, Heidelberg. This work has made use of data from the European Space Agency (ESA) mission Gaia (http://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, http://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. Part of this research was carried out at the Jet Propulsion Laboratory, Caltech, under a contract with NASA. BGM simulations were executed on computers from the Utinam Institute of the Université de Franche-Comté supported by the Région de Franche-Comté and Institut des Sciences de l'Univers (INSU).

J.G. designed BANYAN Σ, compiled the bona fide members, wrote the IDL codes, wrote the manuscript, generated figures and led all analyses; E.E.M. shared parts of the young association literature data, characteristics and bona fide members, and provided general comments; O.R.L. wrote the initial Python translation of BANYAN Σ; L.M., A.R., and D.R. performed the BANYAN I, LACEwING and convergent point tool calculations used in Section 8, respectively; A.R. also provided help with the bona fide members compilation; D.L. and J.K.F. shared ideas and comments; L.P. shared ideas and provided comments especially for the ROC curves analysis; A.R. performed custom Besançon Galactic simulations and wrote the second pagraph of Section 6; and R.D. shared comments and supervised O.R.L.

**Software:** BANYAN Σ (this paper), LACEwING (Riedel et al. 2017b), BANYAN II (Gagné et al. 2014), BANYAN I (Malo et al. 2013), the convergent point tool (Rodríguez et al. 2011), Notability by Ginger Labs, Sublime Text.

## Appendix A

### Coordinate Transformation of the Bayesian Likelihood

In this section, a change of coordinates is applied to the Bayesian likelihood $P_o(\{O_i\}|H_k)$ from the direct observables frame of reference $\{O_i\}$ (sky position, proper motion, etc.) to the Galactic position (XYZ) and space velocity (UVW) reference frame $\{Q_i\}$:

$$P_o(\{O_i\}|H_k) \prod_j dO_j = P_q(\{Q_i\}|H_k) \prod_j dQ_j.$$

This change of coordinates can be expressed in the following form:

$$P_o(\{O_i\}|H_k) = |\tilde{J}| \cdot P_q(\{Q_i\}|H_k), \quad (23)$$

$$\text{with } J_{lm} = \frac{\partial Q_l}{\partial O_m}, \quad (24)$$
where $\tilde{J}$ is a Jacobian matrix, which can be expressed as a
$2 \times 2$ block matrix of four $3 \times 3$ sub-matrices:

$$
\begin{align*}
\tilde{J} &= \begin{bmatrix} \mathcal{A} & \mathcal{B} \\ \mathcal{C} & \mathcal{D} \end{bmatrix}, \\
\mathcal{A} &= \begin{bmatrix} \frac{\partial X}{\partial \alpha} & \frac{\partial X}{\partial \delta} & \frac{\partial X}{\partial \bar{\varpi}} \\ \frac{\partial Y}{\partial \alpha} & \frac{\partial Y}{\partial \delta} & \frac{\partial Y}{\partial \bar{\varpi}} \\ \frac{\partial Z}{\partial \alpha} & \frac{\partial Z}{\partial \delta} & \frac{\partial Z}{\partial \bar{\varpi}} \end{bmatrix}, \\
\mathcal{B} &= \begin{bmatrix} \frac{\partial X}{\partial \mu_{\alpha}} & \frac{\partial X}{\partial \mu_{\delta}} & \frac{\partial X}{\partial \nu} \\ \frac{\partial Y}{\partial \mu_{\alpha}} & \frac{\partial Y}{\partial \mu_{\delta}} & \frac{\partial Y}{\partial \nu} \\ \frac{\partial Z}{\partial \mu_{\alpha}} & \frac{\partial Z}{\partial \mu_{\delta}} & \frac{\partial Z}{\partial \nu} \end{bmatrix}, \\
\mathcal{C} &= \begin{bmatrix} \frac{\partial U}{\partial \alpha} & \frac{\partial U}{\partial \delta} & \frac{\partial U}{\partial \bar{\varpi}} \\ \frac{\partial V}{\partial \alpha} & \frac{\partial V}{\partial \delta} & \frac{\partial V}{\partial \bar{\varpi}} \\ \frac{\partial W}{\partial \alpha} & \frac{\partial W}{\partial \delta} & \frac{\partial W}{\partial \bar{\varpi}} \end{bmatrix}, \\
\mathcal{D} &= \begin{bmatrix} \frac{\partial U}{\partial \mu_{\alpha}} & \frac{\partial U}{\partial \mu_{\delta}} & \frac{\partial U}{\partial \nu} \\ \frac{\partial V}{\partial \mu_{\alpha}} & \frac{\partial V}{\partial \mu_{\delta}} & \frac{\partial V}{\partial \nu} \\ \frac{\partial W}{\partial \mu_{\alpha}} & \frac{\partial W}{\partial \mu_{\delta}} & \frac{\partial W}{\partial \nu} \end{bmatrix}.
\end{align*}
$$

From the definition of Galactic coordinates (Johnson &
Soderblom 1987):

$$
(X, Y, Z) = \varpi \lambda(\alpha, \delta) \quad (25)
$$

it is apparent that **B** = 0. The determinant of **J̄** can be obtained
from the following property of block matrices, and further
simplified:

$$
\begin{align*}
|\bar{J}| &= |\mathcal{A}\mathcal{D} - 0\mathcal{C}|, \\
&= |\mathcal{A}\mathcal{D}|, \\
&= |\mathcal{A}| \cdot |\mathcal{D}|.
\end{align*}
$$

Using the definition of Galactic coordinates in Equation (25),
the sub-matrix **A** can be simplified to:

$$
\begin{align*}
\mathcal{A} &= \begin{bmatrix} \frac{\partial \lambda_0}{\partial \alpha} & \frac{\partial \lambda_0}{\partial \delta} & \lambda_0 \\ \frac{\partial \lambda_1}{\partial \alpha} & \frac{\partial \lambda_1}{\partial \delta} & \lambda_1 \\ \frac{\partial \lambda_2}{\partial \alpha} & \frac{\partial \lambda_2}{\partial \delta} & \lambda_2 \end{bmatrix}, \\
|\mathcal{A}| &= \varpi^2 \lambda_0 \left( \frac{\partial \lambda_1}{\partial \alpha} \frac{\partial \lambda_2}{\partial \delta} - \frac{\partial \lambda_1}{\partial \delta} \frac{\partial \lambda_2}{\partial \alpha} \right) \\
&\quad - \varpi^2 \lambda_1 \left( \frac{\partial \lambda_0}{\partial \alpha} \frac{\partial \lambda_2}{\partial \delta} - \frac{\partial \lambda_0}{\partial \delta} \frac{\partial \lambda_2}{\partial \alpha} \right) \\
&\quad + \varpi^2 \lambda_2 \left( \frac{\partial \lambda_0}{\partial \alpha} \frac{\partial \lambda_1}{\partial \delta} - \frac{\partial \lambda_0}{\partial \delta} \frac{\partial \lambda_1}{\partial \alpha} \right), \\
|\mathcal{A}| &= \varpi^2 f(\alpha, \delta),
\end{align*}
$$

where $f(\alpha, \delta)$ is a function of sky coordinates only.

From the definition of space velocity (see Section 3.2):

$$
(U, V, W) = \varpi N(\alpha, \delta, \mu_{\alpha}, \mu_{\delta}) + \nu M(\alpha, \delta),
$$

where the term cos δ has been omitted in μα cos δ, the sub-
matrix ⃗D can be simplified to:

$$
\begin{align*}
\mathcal{D} &= \begin{bmatrix}
\frac{\varpi}{\partial \mu_{\alpha}} \quad \frac{\varpi}{\partial \mu_{\delta}} & M_0 \\
\frac{\varpi}{\partial \mu_{\alpha}} \quad \frac{\varpi}{\partial \mu_{\delta}} & M_1 \\
\frac{\varpi}{\partial \mu_{\alpha}} \quad \frac{\varpi}{\partial \mu_{\delta}} & M_2
\end{bmatrix}, \\
|\mathcal{D}| &= \varpi^2 M_0 \left( \frac{\partial N_1}{\partial \mu_{\alpha}} \frac{\partial N_2}{\partial \mu_{\delta}} - \frac{\partial N_1}{\partial \mu_{\delta}} \frac{\partial N_2}{\partial \mu_{\alpha}} \right) \\
&\quad - \varpi^2 M_1 \left( \frac{\partial N_0}{\partial \mu_{\alpha}} \frac{\partial N_2}{\partial \mu_{\delta}} - \frac{\partial N_0}{\partial \mu_{\delta}} \frac{\partial N_2}{\partial \mu_{\alpha}} \right) \\
&\quad + \varpi^2 M_2 \left( \frac{\partial N_0}{\partial \mu_{\alpha}} \frac{\partial N_1}{\partial \mu_{\delta}} - \frac{\partial N_0}{\partial \mu_{\delta}} \frac{\partial N_1}{\partial \mu_{\alpha}} \right), \\
|\mathcal{D}| &= \varpi^2 g(\alpha, \delta).
\end{align*}
$$

where g(α, δ) is a function of sky coordinates only. The fact that
g does not depend on the proper motion components arises from
the fact that all components of **N**, defined in Section 3.2, depend
linearly on the proper motion components. It follows that:

$|
\bar{J}
| \propto \varpi^4
.$

**Appendix B**
Solving the Marginalization Integrals

In this section, the analytical solution to the marginalization
integrals of Equation (3), over distance and radial velocity, is
developed. The index *k* that referring to hypothesis *H*<sub>*k*</sub> is
ignored in this section for simplicity. The Bayesian likelihood
in the Galactic frame of reference {*Q*<sub>*i*</sub>} described in
Equation (11) can be inserted in Equation (3) to obtain:

$$
\mathcal{P}(\{O_i\} | H) = \int_0^\infty \int_{-\infty}^\infty \varpi^4 \frac{e^{-\frac{1}{2}M^2}}{\sqrt{(2\pi)^6 |\bar{\Sigma}|}} d\nu \, d\bar{\varpi},
$$

and the Mahalanobis distance $\mathcal{M}$ defined in Equation (10) can
be used to develop it further:

$$
\mathcal{P}(\{O_i\} | H) = C_0 \int_0^\infty I_0(\varpi) \varpi^4 e^{-\frac{1}{2}(\bar{\Gamma}, \bar{\Gamma})\varpi^2 + (\bar{\Gamma}, \bar{\tau})\varpi} d\varpi, \quad (26)
$$

$$
\begin{align}
I_0(\varpi) &= \int_{-\infty}^{\infty} e^{-\frac{1}{2} \langle \bar{\Omega}, \bar{\Omega} \rangle \nu^2 + \langle \bar{\Omega}, \bar{\tau} \rangle - \langle \bar{\Omega}, \bar{\Gamma} \rangle \nu} d\nu, \nonumber \\
C_0 &= \frac{e^{-\frac{1}{2} \langle \bar{\tau}, \bar{\tau} \rangle}}{\sqrt{(2\pi)^6 |\bar{\Sigma}|}}. \tag{27}
\end{align}
$$

Equation (27) can be solved with the identity:

$$
\int_{-\infty}^{\infty} e^{-ax^2 - bx} dx = \sqrt{\frac{\pi}{a}} e^{b^2/4a},
$$
with $a = \frac{1}{2} \langle \bar{\Omega}, \bar{\Omega} \rangle$ and $b = - \langle \bar{\Omega}, \bar{\tau} \rangle + \langle \bar{\Omega}, \bar{\Gamma} \rangle \varpi$.
The term in $b^2$ can be developed into a second-degree polynomial in $\varpi$:

$$ \frac{b^2}{4a} = \frac{1}{2} \frac{\langle \bar{\Omega}, \bar{\Gamma} \rangle^2}{\langle \bar{\Omega}, \bar{\Omega} \rangle} \varpi^2 - \frac{\langle \bar{\Omega}, \bar{\Gamma} \rangle \langle \bar{\Omega}, \bar{\tau} \rangle}{\langle \bar{\Omega}, \bar{\Omega} \rangle} \varpi + \frac{1}{2} \frac{\langle \bar{\Omega}, \bar{\tau} \rangle^2}{\langle \bar{\Omega}, \bar{\Omega} \rangle}, $$

which can be inserted back into Equation (26):

$$ P(\{O_i\}|H) = \frac{|\bar{\Omega}|^{-1}e^{-\zeta}}{\sqrt{(2\pi)^5|\bar{\Sigma}|}} \int_0^{\infty} \varpi^4 e^{-\beta\varpi^2 - \gamma\varpi} d\varpi, $$

$$ \begin{aligned} \text{where } \beta &= \frac{\langle \bar{\Gamma}, \bar{\Gamma} \rangle}{2} - \frac{1}{2} \frac{\langle \bar{\Omega}, \bar{\Gamma} \rangle^2}{\langle \bar{\Omega}, \bar{\Omega} \rangle}, \\ \gamma &= \frac{\langle \bar{\Omega}, \bar{\Gamma} \rangle \langle \bar{\Omega}, \bar{\tau} \rangle}{\langle \bar{\Omega}, \bar{\Omega} \rangle} - \langle \bar{\Gamma}, \bar{\tau} \rangle, \\ \zeta &= \frac{\langle \bar{\tau}, \bar{\tau} \rangle}{2} - \frac{1}{2} \frac{\langle \bar{\Omega}, \bar{\tau} \rangle^2}{\langle \bar{\Omega}, \bar{\Omega} \rangle}. \end{aligned} \quad (28) $$

The solution to this integral is given by Erdelyi (1955) and Gradshteyn & Ryzhik (2014):

$$ \int_0^\infty x^n e^{-\beta x^2 - \gamma x} dx = \frac{n! e^{\gamma^2/8\beta}}{(2\beta)^{(n+1)/2}} \mathcal{D}_{-(n+1)}(\gamma/\sqrt{2\beta}), $$

where $\mathcal{D}_m(x)$ is a parabolic cylinder function (Magnus & Oberhettinger 1948).

The $m = -5$ ($n = 4$) corresponds to the required integral, and the corresponding parabolic cylinder function can be developed as:

$$ \mathcal{D}_{-5}(x) = \frac{e^{x^2/4}}{24} \left( \sqrt{\frac{\pi}{2}} (x^4 + 6x^2 + 3) \operatorname{erfc} \left( \frac{x}{\sqrt{2}} \right) - (x^3 + 5x)e^{-x^2/2} \right). $$

Equation (28) becomes:

$$ \mathcal{P}(\{O_i\}|H) = \frac{3 \mathcal{D}_{-5}(\gamma/\sqrt{2\beta}) e^{\gamma^2/8\beta - \zeta}}{4 |\bar{\Omega}| \sqrt{\pi^5 \beta^5 |\bar{\Sigma}|}}. \quad (29) $$

The term in $e^{x^2/4}$ in the definition of $\mathcal{D}_{-5}(x)$ can become very large for typical values of $x$, making the numerical computation of $\mathcal{D}_{-5}(x)$ unstable. To avoid this problem, a modified parabolic cylinder function $\mathcal{D}'_{-5}(x)$ can be defined so that the large term is combined with the exponential term in Equation (29):

$$ \begin{aligned} \mathcal{P}(\{O_i\}|H) &= \frac{1}{32} \frac{\mathcal{D}'_{-5}(\gamma/\sqrt{2\beta}) e^{\gamma^2/4\beta - \zeta}}{|\bar{\Omega}|\sqrt{\pi^5 \beta^5 |\bar{\Sigma}|}}, \\ \mathcal{D}'_{-5}(x) &= 24 e^{-x^2/4} \mathcal{D}_{-5}(x). \end{aligned} $$

which completes the analytical solution of the Bayesian likelihood. The $1/32$ factor will be ignored here because it will disappear in the marginalization of the Bayesian likelihood (Equation (1)), and other multiplicative factors independent of the Bayesian hypothesis have already been ignored in determining the Jacobian of the coordinate transformation (see Appendix A).

## Appendix C
### Determining the Optimal Radial Velocity and Distance

The optimal radial velocity $\nu_o$ and distance $\varpi_o$ that maximize the value of the non-marginalized Bayesian likelihood $\mathcal{P}_o(\{O_i\}|H)$ can be obtained by solving the system of

equations:

$$ \begin{aligned} \frac{\partial \ln \mathcal{P}_o(\{O_i\}|H)}{\partial \nu} \Bigg|_{\nu=\nu_o, \varpi=\varpi_o} &= 0, \\ \frac{\partial \ln \mathcal{P}_o(\{O_i\}|H)}{\partial \varpi} \Bigg|_{\nu=\nu_o, \varpi=\varpi_o} &= 0, \end{aligned} $$

which can be developed with Equation (10):

$$ \begin{aligned} 0 &= \langle \bar{\Omega}, \bar{\Omega} \rangle \nu_o + \langle \bar{\Omega}, \bar{\Gamma} \rangle \varpi_o - \langle \bar{\Omega}, \bar{\tau} \rangle, \\ 0 &= \langle \bar{\Gamma}, \bar{\Gamma} \rangle \varpi_o^2 + \langle \bar{\Omega}, \bar{\Gamma} \rangle \varpi_o \nu_o - \langle \bar{\Gamma}, \bar{\tau} \rangle \varpi_o - 4. \end{aligned} $$

This system of equations has two solutions:

$$ \begin{aligned} \varpi_o &= \frac{-\gamma \pm \sqrt{\gamma^2 + 32\beta}}{4\beta}, \\ \nu_o &= \frac{4 + \langle \bar{\Gamma}, \bar{\tau} \rangle \varpi_o - \langle \bar{\Gamma}, \bar{\Gamma} \rangle \varpi_o^2}{\langle \bar{\Omega}, \bar{\Gamma} \rangle \varpi_o}, \end{aligned} $$

where $\gamma$ and $\beta$ are defined in Equations (14) and (13).

Any combination of $\gamma$ and $\beta$ that respects the following inequality:

$$ \begin{aligned} \sqrt{1 + 32\beta/\gamma^2} &> 1, \\ \text{i.e., } &\beta > 0, \end{aligned} $$

will yield an unphysical negative distance for the negative root of $\varpi_o$. Since multivariate Gaussians have $\beta > 0$ by definition, the inequality is always respected. As a consequence, only the positive root of $\varpi_o$ has a physical meaning.

Error bars on the optimal radial velocity $\sigma_\nu$ and distance $\sigma_\varpi$ can be defined by measuring the characteristic width of $\mathcal{P}_o(\{O_i\}|H)$ along $\nu$ and $\varpi$ in the vicinity of $(\nu_o, \varpi_o)$. The effect of the Jacobian term $\varpi^4$ will be ignored to obtain an analytical approximation of $(\sigma_\nu, \sigma_\varpi)$.

The relation between the expectancy $E(x)$ of a variable and the characteristic width of a Gaussian function $G(x)$ can be used to determine $\sigma_\nu$ and $\sigma_\varpi$:

$$ \begin{aligned} \sigma_x &= \sqrt{E(x^2) - E(x)^2}, \\ E(x) &= \int_{-\infty}^{\infty} xG(x)dx. \end{aligned} $$

In the case of $\sigma_\nu$, this yields:

$$ \begin{aligned} E(\nu) &= \frac{\int_{-\infty}^{\infty} \nu e^{-\beta_{\nu} \nu^2 - \gamma_{\nu} \nu} d\nu}{\int_{-\infty}^{\infty} e^{-\beta_{\nu} \nu^2 - \gamma_{\nu} \nu} d\nu}, \\ &= \frac{\gamma_{\nu}}{2\beta_{\nu}}, \\ E(\nu^2) &= \frac{\int_{-\infty}^{\infty} \nu^2 e^{-\beta_{\nu} \nu^2 - \gamma_{\nu} \nu} d\nu}{\int_{-\infty}^{\infty} e^{-\beta_{\nu} \nu^2 - \gamma_{\nu} \nu} d\nu}, \\ &= \frac{2\beta_{\nu} + \gamma_{\nu}^2}{4\beta_{\nu}^2}, \\ \beta_{\nu} &= \frac{\langle \bar{\Omega}, \bar{\Omega} \rangle}{2}, \\ \gamma_{\nu} &= \varpi_o \langle \bar{\Omega}, \bar{\Gamma} \rangle - \langle \bar{\Omega}, \bar{\tau} \rangle, \end{aligned} $$
leading to:

$$\sigma_{\nu} = \frac{1}{\sqrt{2\beta_{\nu}}},$$

$$\sigma_{\nu} = |\hat{\Omega}|^{-1}.$$

The case of $\sigma_{\bar{\omega}}$ requires the introduction of a new variable $\bar{\omega}'$ that is defined in the range $]-\infty, \infty[$ and matches the distance $\bar{\omega}' = \bar{\omega}$ for $\bar{\omega}' \ge 0$. Assuming that $\bar{\omega}_0 \gg \sigma_{\bar{\omega}}$ will ensure that the Bayesian likelihood $P_0(\nu, \bar{\omega}') \approx 0$ for all negative values of $\bar{\omega}'$. It follows that:

$$E(\bar{\omega}) \approx \frac{\int_{-\infty}^{\infty} \bar{\omega}' e^{-\beta_{\bar{\omega}}\bar{\omega}^2 - \gamma_{\bar{\omega}}\bar{\omega}} d\bar{\omega}'}{\int_{-\infty}^{\infty} e^{-\beta_{\bar{\omega}}\bar{\omega}^2 - \gamma_{\bar{\omega}}\bar{\omega}} d\bar{\omega}'} \\ = \frac{\gamma_{\bar{\omega}}}{2\beta_{\bar{\omega}}}$$

$$E(\bar{\omega}^2) \approx \frac{2\beta_{\bar{\omega}} + \gamma_{\bar{\omega}}^2}{4\beta_{\bar{\omega}}^2},$$

$$\beta_{\bar{\omega}} = \frac{\langle \bar{\Gamma}, \bar{\Gamma} \rangle}{2},$$

$$\gamma_{\bar{\omega}} = \nu_0 \langle \hat{\Omega}, \bar{\Gamma} \rangle - \langle \bar{\Gamma}, \bar{\tau} \rangle,$$

leading to:

$$\sigma_{\bar{\omega}} \approx |\bar{\Gamma}|^{-1}.$$

The optimal distance and radial velocity do not correspond exactly to the statistical distance and radial velocities defined in the BANYAN II formalism (Gagné et al. 2014). The latter are obtained by maximizing the Bayesian likelihood in one dimension after the other dimension was marginalized. The optimal distance and radial velocity maximize the Bayesian probability of a given hypothesis as a couple, whereas the BANYAN II statistical distance maximizes the Bayesian probability when radial velocity is treated as an unknown parameter, and vice versa.

## ORCID iDs

Jonathan Gagné https://orcid.org/0000-0002-2592-9612

Eric E. Mamajek https://orcid.org/0000-0003-2008-1488

Lison Malo https://orcid.org/0000-0002-8786-8499

Adric Riedel https://orcid.org/0000-0003-1645-8596

David Rodriguez https://orcid.org/0000-0003-1286-5231

David Lafrenière https://orcid.org/0000-0002-6780-4252

Jacqueline K. Faherty https://orcid.org/0000-0001-6251-0573

Olivier Roy-Loubier https://orcid.org/0000-0003-2005-5626

Laurent Pueyo https://orcid.org/0000-0003-3818-408X

Annie C. Robin https://orcid.org/0000-0001-8654-9499

René Doyon https://orcid.org/0000-0001-5485-4675

## References

* Abt, H. A. 1981, ApJS, 45, 437

* Abt, H. A. 2004, ApJS, 155, 175

* Abt, H. A. 2008, ApJS, 176, 216

* Abt, H. A. 2009, ApJS, 180, 117

* Abt, H. A., & Levato, H. 1978, PASP, 90, 201

* Abt, H. A., & Levy, S. G. 1985, ApJS, 59, 229

* Abt, H. A., & Morrell, N. I. 1995, ApJS, 99, 135

* Adams, W. S., Joy, A. H., Humason, M. L., & Brayton, A. M. 1935, ApJ, 81, 187

* Alam, S., Albareti, F. D., Allende Prieto, C., et al. 2015, ApJS, 219, 12

* Alcalá, J. M., Covino, E., Torres, G., et al. 2000, A&A, 353, 186

* Aller, K. M., Liu, M. C., Magnier, E. A., et al. 2016, ApJ, 821, 120

* Allers, K. N., Gallimore, J. F., Liu, M. C., & Dupuy, T. J. 2016, ApJ, 819, 133

* Allers, K. N., & Liu, M. C. 2013, ApJ, 772, 79

* Allison, R. J., Goodwin, S. P., Parker, R. J., et al. 2009, MNRAS, 395, 1449

* Alves de Oliveira, C., Moraux, E., Bouvier, J., & Bouy, H. 2012, A&A, 539, 151

* Anderson, E., & Francis, C. 2012, AstL, 38, 331

* Ansdel, M., Gaidos, E., Rappaport, S. A., et al. 2016, ApJ, 816, 69

* Appenzeller, I., Reitermann, A., & Stahl, O. 1988, PASP, 100, 815

* Arenou, F. 2011, in AIP Conf. Proc. 1346, International Workshop on Double and Multiple Stars: Dynamics, ed. J. A. Docobo, V. S. Tamazian, & Y. Y. Balega (Melville, NY: AIP), 107

* Barrado y Navascués, D., Stauffer, J. R., & Jayawardhana, R. 2004, ApJ, 614, 386

* Beavers, W. I., & Cook, D. B. 1980, ApJS, 44, 489

* Bell, C. P. M., Mamajek, E. E., & Naylor, T. 2015, MNRAS, 454, 593

* Bell, C. P. M., Murphy, S. J., & Mamajek, E. E. 2017, MNRAS, 468, 1198

* Benedict, G. F., Tanner, A. M., Cargile, P. A., & Ciardi, D. R. 2014, AJ, 148, 108

* Bergbusch, P. A., & Vandenberg, D. A. 1992, ApJS, 81, 163

* Bertelli, G., Bressan, A., Chiosi, C., Fagotto, F., & Nasi, E. 1994, A&AS, 106, 275

* Bertelli, G., Girardi, L., Marigo, P., & Nasi, E. 2008, A&A, 484, 815

* Bertelli, G., Nasi, E., Girardi, L., & Marigo, P. 2009, A&A, 508, 355

* Bilíková, J., Chu, Y.-H., Gruendl, R. A., & Maddox, L. A. 2010, AJ, 140, 1433

* Binnendijk, L. 1946, AnLei, 19, B1

* Bishop, C. M., & Nasrabadi, N. M. 2007, JEI, 16, 049901

* Blaauw, A. 1946, PGro, 52, 1

* Blake, C. H., Charbonneau, D., & White, R. J. 2010, ApJ, 723, 684

* Bobylev, V. V., & Bajkova, A. T. 2007, AstL, 33, 571

* Bobylev, V. V., Goncharov, G. A., & Bajkova, A. T. 2006, ARep, 50, 733

* Bonnefoy, M., Boccaletti, A., Lagrange, A.-M., et al. 2013, A&A, 555, A107

* Bourgés, L., Lafrasse, S., Mella, G., et al. 2014, in ASP Conf. Ser. 485, Astronomical Data Analysis Software and Systems XXIII, ed. N. Manset & P. Forshay (San Francisco, CA: ASP), 223

* Bouvier, J., & Appenzeller, I. 1992, A&AS, 92, 481

* Bouy, H., & Martín, E. L. 2009, A&A, 504, 981

* Bowler, B. P., Liu, M. C., Mawet, D., et al. 2017, AJ, 153, 18

* Bragança, G. A., Daflon, S., Cunha, K., et al. 2012, AJ, 144, 130

* Brandt, T. D., & Huang, C. X. 2015, ApJ, 807, 24

* Brandt, T. D., Kuzuhara, M., McElwain, M. W., et al. 2014, ApJ, 786, 1

* Breger, M. 1984, A&AS, 57, 217

* Burgasser, A. J., Lopez, M. A., Mamajek, E. E., et al. 2016, ApJ, 820, 32

* Buscombe, W. 1969, MNRAS, 144, 31

* Bystrov, N. F., Polojentsev, D. D., Potter, H. I., et al. 1994, BICDS, 44, 3

* Carmona, A., van den Ancker, M. E., & Henning, T. 2007, A&A, 464, 687

* Carpenter, J. M., Mamajek, E. E., Hillenbrand, L. A., & Meyer, M. R. 2006, ApJ, 651, L49

* Cartwright, A., & Whitworth, A. P. 2004, MNRAS, 348, 589

* Casewell, S. L., Jameson, R. F., & Dobbie, P. D. 2006, MNRAS, 365, 447

* Casewell, S. L., Littlefair, S. P., Burleigh, M. R., & Roy, M. 2014, MNRAS, 441, 2644

* Cayrel de Strobel, G., Soubiran, C., & Ralite, N. 2001, A&A, 373, 159

* Cenarro, A. J., Cardiel, N., Vazdekis, A., & Gorgas, J. 2009, MNRAS, 396, 1895

* Cenarro, A. J., Peletier, R. F., Sánchez-Blázquez, P., et al. 2007, MNRAS, 374, 664

* Chabrier, G., & Baraffe, I. 1997, A&A, 327, 1039

* Chen, C. H., Mamajek, E. E., Bitner, M. A., et al. 2011, ApJ, 738, 122

* Christy, J. W., & Walker, R. L. J. 1969, PASP, 81, 643

* Chubak, C., & Marcy, G. 2011, in AAS Meeting, 217, 434.12

* Chubak, C., Marcy, G., Fischer, D. A., et al. 2012, arXiv:1207.6212

* Cieza, L., Padgett, D. L., Stapelfeldt, K. R., et al. 2007, ApJ, 667, 308

* Cohen, M., & Kuhi, L. V. 1979, ApJS, 41, 743

* Corbally, C. J. 1984, ApJS, 55, 657

* Corporon, P., Lagrange, A.-M., & Beust, H. 1996, A&A, 310, 228

* Cottaar, M., Covey, K. R., Foster, J. B., et al. 2015, ApJ, 807, 27

* Covino, E., Alcalá, J. M., Allain, S., et al. 1997, A&A, 328, 187

* Cowley, A., Cowley, C., Jaschek, M., & Jaschek, C. 1969, AJ, 74, 375

* Cowley, A., & Fraquelli, D. 1974, PASP, 86, 70

* Cruz, K. L., Reid, I. N., Kirkpatrick, J. D., et al. 2007, AJ, 133, 439

* Cucchiaro, A., Jaschek, M., Jaschek, C., & Macau-Hercot, D. 1976, A&AS, 26, 241