BANYAN. XI. The BANYAN $\Sigma$ multivariate Bayesian algorithm to identify members of young associations within 150 pc

BANYAN $\Sigma$ 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, modelled with multivariate Gaussians in 6-dimensional 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 $\rho$ 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\c{c}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 $\Sigma$ benefits from an analytical solution to the Bayesian marginalization integrals 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 $\Sigma$ achieves a better classification performance than other moving group classification tools, 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, 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 $\Sigma$ are made available with this publication. (shortened)

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 co-eval populations, making them valuable to measure the initial mass function and serve as age calibrators for members across all masses.
Recent surveys (e.g., Gagné et al. 2014Gagné et al. , 2015bKellogg et al. 2015;Gagné et al. 2015c;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. 2001aZuckerman et al. , 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 U V W , with typical velocity dispersions below ∼ 3 km s −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 motion are available. These include the convergent point tool (e.g., Mamajek 2005;Torres et al. 2006;Rodriguez et al. 2011), various goodnessof-fit metrics (e.g. Kraus et al. 2014;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 AssociatioNs), an algorithm based on Bayesian infer-ence where moving groups are modelled with unidimensional Gaussian distributions in Galactic coordinates XY Z and space velocities U V W . 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 modelled 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 W 1 and H − W 2 versus M W 1 ) 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. 2013Malo et al. , 2014Gagné et al. 2015c;Faherty et al. 2016;Liu et al. 2016). The majority of these classification tools include only the seven youngest  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 Σ, 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 welldefined associations within ∼ 150 pc 1 . BANYAN Σ in-cludes 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 multivariate Gaussians that capture full correlations in the XY ZU V W distribution of members, including those in mixed spatial-kinematic coordinates. Two versions of the BANYAN Σ code (IDL 2 and python 3 ) are made publicly available (Gagné et al. 2018b,a), and a web portal is made available for singleobject queries 4 .
Most previously available classification tools rely on time-consuming algorithms, such as numerical integrals, which make it challenging to analyze large datasets 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 Σ, 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 Σ is also the first classification tool to include the Taurus, ρ Ophiuchi and Corona Australis star-forming regions (e.g., Wichmann et al. 2000;Reipurth 2008), 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 latter, 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 Σ, 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 Σ in Sec-tion 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 Σ 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 Σ. This work is concluded in Section 10.

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 Σ 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 (α,δ), proper motion (µ α ,µ δ ), radial velocity (ν) and distance ( ), 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: 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 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: It is often the case that the radial velocity (ν) and/or distance ( ) 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: 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 XY Z and space velocity U V W 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: where {O i } represents the set of observables excluding the radial velocity ν and distance . Equation (4) inherently ignores the Jacobian of the transformation {O i } → {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 ν and , covering −35 to 35 km s −1 and 0.1 to 200 pc, respectively. On each point (ν, ) 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 } , ν, )}|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 (|ν| > 35 km s −1 ) or distant ( > 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 XY ZU V W 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 Σ. 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: where P m (µ α , µ δ ) 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 XY Z and U V W space, which are freely rotated along any axis. Models for 7 young associations were included: TW Hya (TWA; de La Reza et al. 1989;Kastner et al. 1997), β Pictoris (β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.

BANYAN Σ: AN IMPROVED ALGORITHM
The framework of BANYAN Σ 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 Σ is described in this section, and additional improvements with respect to the kinematic models are described in Sections 4, 5 and 6. The models of BANYAN Σ 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 Σ is analogous to a Bayesian classification algorithm with a Gaussian mixture model (e.g., see Bishop &Nasrabadi 2007 andMcLachlan &Peel 2000).
The multivariate Gaussian models of BANYAN Σ 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 ν and distance 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 Σ, 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 to 3.9 describe additional options to the BANYAN Σ algorithm, to include measurements of radial velocity and/or distance, constraints from spectrophotometric observables, or to ignore the spatial distribution of young associations.

Kinematic Models
The distribution of stars in young stellar associations and the Galactic field are modelled with multivariate Gaussian distributions in XY ZU V W space. This is a generalization over the freely rotating individual XY Z and U V W 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: is the Mahalanobis distance and the bar (e.g.,v) and double-bar (e.g.,M ) symbols are used to indicate vectors and matrices, respectively, in XY ZU V W space.Q is a 6-dimensional vector built from {Q i }, which correspond to the XY ZU V W 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τ vector and indicate the center of the association; six are stored in the diagonal of the covariance matrixΣ and indicate the 6D size of the association; and 15 more are stored in the independent off-diagonal elements ofΣ and indicate the orientation of the ellipsoid in 6D space, or equivalently the correlations between each combination of coordinates.x T indicates a vector transposition and s −1 a matrix inverse.
The Bayesian likelihoods in the Galactic frame of reference {Q i } can thus be written as: where the dependencies ofτ andΣ on the association index k are implicit.
A simple approach for obtaining the parameters of a kinematic model is to calculate the average positionτ of the members in XY ZU V W space and their variances and covariances to build theΣ matrix. A method that is more robust to outliers, and accounts for individual measurement errors, is presented in Section 5.

Change of Coordinates
Solving the Bayesian likelihood in Equation (3) requires applying a change of coordinates from the observables frame of reference {O i } to the Galactic frame of reference {Q i }. The equations for this transformation are detailed by Johnson & Soderblom (1987), where the components ofQ i can be written as a linear combination of the radial velocity ν and the distance : the components ofΩ andΓ are: and symbols in bold represent 3D vectors or matrices in XY Z or U V W space. The vectors λ, M and N transform the sky position, proper motion, radial velocity and distance to XY ZU V W following: where l and b are the Galactic longitude and latitude, α and δ are the right ascension and declination, µ α cos δ and µ δ are the proper motion, and κ ≈ 4.74 · 10 −3 corresponds to 10 −3 AU/yr so that proper motions are expressed in mas yr −1 . The 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 T where the first row has the opposite sign of that defined by Johnson & Soderblom (1987), so that U points towards the galactic center and U V W forms a right-handed system 5 : Equation (6) can be used to express 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 P q → P o . This step is detailed in Appendix A, and yields: Inserting the coordinate transformation (Equation 6) into the kinematic model defined in Equation (5) yields: All terms in Ω i , Γ i and τ i,k can be described as scalar products induced by the inverse covariance matrix σ −1 ij,k , such as: 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 simplifications we get: This completes the coordinate change of the Bayesian likelihood: 5 This display of the T matrix follows the column major convention (e.g., IDL or Forftran). Because Python follows the row major matrix convention, T must be transposed in that language.

Solving the Marginalization Integrals
A complete analytical solution to Equation (3) is developed in Appendix B and yields: 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Ω,Γ,τ (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 6 components of the 6D vectors, as well as β, γ, ζ and P({O i }|H) for the full array of stars at once.

Optimal Radial Velocity and Distance
The optimal values for the radial velocity and distance (ν o , 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: where σ and σ ν represent statistical 1σ error bars on the optimal values. These two values are given for each hypothesis and each star as an output of BANYAN Σ, 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 (ν o , 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 σ away from the prediction. The values of σ and σ ν will therefore be unreliable in such a situation.

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 ν o and distance 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 XY Z and space velocity U V W in the vicinity of the maximum of the membership probability distribution function. The sky position, proper motion, ν o , o and proper motion errors are propagated in XY ZU V W space to obtain a Galactic error vectorσ q = (σ X , σ Y , σ Z , σ U , σ V , σ W ).
It is then possible to add these errors in quadrature to the diagonal elements ofΣ without affecting the orientation of the multivariate Gaussian with: where diag − M extracts the diagonal elements of ma-trixM and diag + (v) builds a diagonal matrix with vectorv.
The need for a time-consuming matrix inversion ofΣ for each star to evaluate Equation (10), can be avoided with:Σ Including an approximated effect of the proper motion measurement errors will thus require a calculation of (1) theΩ,Γ andτ vectors and their scalar products, (2) ν o and o , (3) theσ q vector and the corresponding inflated matrixΣ , and (4) all quantities from the beginning of Section 3.2 obtained with the updated scalar product based onΣ −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 XY ZU V W 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.

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 Σ formalism, as long as the parallax factors (ψ α , ψ δ ), 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 (∆ απ , ∆ δπ ) of a star is given by Smart & Green (1977, Chap. 9, p. 221): ∆ δπ = cos δ sin e sin s − cos α sin δ cos s − sin α sin δ cos e sin s , where e is the obliquity of the ecliptic of the Earth's orbit and s is the ecliptic longitude of the Sun at a given epoch 6 . These equations can be simplified by grouping all epoch-dependent terms into (φ α , φ δ ): where t is the epoch. The apparent motion (µ α , µ δ ) of a star between epochs t 1 and t 2 will thus be given by: = µ δ + ψ δ , and similarly : Since Equation (6) is linear in proper motion components, it can be expressed as a function of apparent motion with the form: whereΓ is a function of apparent motion (µ α cos δ, µ δ ), the quantity that is directly measured, instead of true proper motion (µ α cos δ, µ δ ). BecauseQ only appears in the Bayesian likelihood as relative to the center of the moving group model τ , the effect of parallax motion can be fully accounted for by shifting the U V W center of the young association kinematic model by +T A ψ, which is equivalent to shiftinḡ τ by +Φ:τ =τ +Φ.
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τ →τ in the BANYAN Σ formalism. This requires measurements of the parallax factors (ψ α , ψ δ ) for each star in addition to measurements of their apparent motion. In practice, this correction is applied by the BANYAN Σ software only when the keyword use psi is explicitly used. This indicates that: (1) the proper motion that is input to BANYAN Σ 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 Σ 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.

Additional Kinematic Observables
Only sky position and proper motion are required for BANYAN Σ 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 ) and BANYAN II (Gagné et al. 2014). In cases where a radial velocity measurement is available, Equation (3) can be rewritten as: where P m (ν) is the probability density function that represents the radial velocity measurement. Assuming that P m (ν) is a Gaussian distribution centered on ν m with a characteristic width of σ ν,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: The case where a distance measurement is available can be solved in a similar way, by replacing: 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 τ ,τ must be cumulated).

Photometric Observables
It is possible to constrain the distance of a star from its position in a color-magnitude or spectral typemagnitude 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 Σ 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 m and its uncertainty σ ,m must be used, one for each hypothesis, because they derive from different colormagnitude 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 ( m , σ ,m ) for the field and each young association, and provide them to BANYAN Σ for a full inclusion of these constrainst in the Bayesian probabilities. Multiple color-magnitude diagrams can also be combined into single measurements of ( m , σ ,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 σ ,m . Such unrealistically precise constraints on the distance of a star would hinder the ability of BANYAN Σ to correctly identify the candidate members of a young association. Both the IDL and Python implementations of BANYAN Σ 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 Σ 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 Σ through the Bayesian priors: the vertical distances in a colormagnitude 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 Σ 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α, 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 Σ prior probabilities. This framework allows users to add observables in the BANYAN Σ membership determination without needing to change its algorithm, and remains accurate as long as no kinematic measure-ments 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 Σ with the keyword ln priors. These data are passed to BANYAN Σ 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 documentations, 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 Σ, but they will be provided in future work as they are developed to target specific types of members.

Ignoring the Galactic Position XY Z
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 XY Z as well as U V W prevents an exploration of that potentially missing population of members.
The BANYAN Σ framework can be adapted to rely on U V W only, by artificially setting the first three diagonal elements of the covariance matricesΣ 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 Σ formalism, and the general solution presented in Equation (12) remains unchanged. An option is provided in BANYAN Σ to ignore the spatial XY Z 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.

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 XY ZU V W space 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 9 associations with incomplete kinematics are described in Section 4.2 and their adopted average distances and radial velocities are listed in Table 4. A description the Argus association, included in BANYAN II but excluded here, is provided in Section 4.3. A few individual ob-jects 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 crossmatched 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;Torres et al. 2006;Bobylev 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).

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. (2017a) 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 comoving 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. (2017b) used a preliminary version of BANYAN Σ to show that the variable T2.5 brown dwarf SIMP J013656.5+093347 is likely a ∼ 13 M Jup 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 co-eval young stars, and a stream of stars with heterogeneous compositions and ages. Soderblom & Mayor (1993) 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) 7 -are astrometric radial velocities assuming moving group membership and an average U V W velocity. These measurements were rejected from our compilation.
The Octans association (OCT; Torres et al. 2008) is a group of young stars at ≈ 120 pc from the Sun, which has not been characterised as well as other young associations mainly due to its sky position located far in the Southern hemisphere (declinations 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 colormagnitude 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.

Associations With Partial Kinematics
This section describes the 9 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. Thes 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 8 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 4.
Instead of assigning the exact same average distance or radial velocity to each members with a missing observable, artifical 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 4). This avoids artifically placing several members along lines in XY Z and U V W spaceor along planes in U V W 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.

Association
Rejection criterion Table 5 continued  (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 χ 1 For). In order to identify more members, the Dias et al. (2014) list of 4 102 XFOR candidates was cross-matched with the Gaia-DR1. Of the 261 matches with a parallax measurement, only 9 are located within 10 pc of the χ 1 For system in XY Z 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 −1 of the star χ 1 For in U V W space if we assume the same radial velocity measurement, and they are therefore included as highlikelihood candidate members. Two additional XFOR members were identified by Alessi et al. (priv. comm.): HD 21434 and HD 17864. Both 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 4 of its members benefit from full kinematic measurements (a Car, HD 76230, H Vel, OY Vel).
ρ Ophiuchi (ROPH) is the nearest star-forming cloud complex to the Sun. It has been the subject of exten-sive studies in the past decades (e.g. see Reipurth et al. 1991;Wilking et al. 2008;Reipurth 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±3 pc. They also estimate an average radial velocity of −6.3±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 Σ 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  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 Σ, 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. Esplin et al. (2014) were adopted here, without differenciating 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 Table 6. New bona fide members with full kinematics compiled in this work.

Designation
Ref a Reference that designated this object as a candidate member of the young association. b 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.
Note-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 detail.
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 2 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) 9 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 ). 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.

Rejected Associations
The Argus association was removed entirely from the models of BANYAN Σ, 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  and Hercules-Lyra associations (Gaidos 1998

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 Σ kinematic models and are listed in Table 10.
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 ± 0.7 km s −1 from the core of CAR and at 5.5 ± 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 βPMG, COL and AB-DMG, because of contradictory measurements for its radial velocity (Kharchenko et al. 2007;Montes et al. 2001;López-Santiago et al. 2010), and Gagné et al. (2014) incorrectly listed it as a bona fide member of both β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 , and Malo et al. (2013) defined it as a bona fide member of COL because it lies nearer in XY Z space, although they note that its membership is ambiguous. Here HD 23524 is excluded from the BANYAN Σ 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 Σ kinematic models.
2MASS J06085283-2753583 was identified by Rice et al. (2010) as a candidate member of β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 proper motion. This object was therefore not included in the kinematic models of BANYAN Σ.

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 6.
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 (5 in βPMG, 8 in THA, 2 in ABDMG and one in COL). The U V W 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 lowmass 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 U V W position that is located at 8.9 km s −1 from the central position of βPMG, and only 3.4 km s −1 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 Riedel et al. 2014); 5/6 were already included in the list of bona fide members presented here. The remaining star, AF Psc, has a parallax measurement by van , 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 to the list of bona fide members.

Calculation of the 6D Kinematics
The Galactic positions XY Z and space velocities U V W , 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 4element Monte Carlo approach was used to propagate error bars in XY ZU V W space, by adopting the standard deviation of each coordinate as its measurement error, and therefore assuming that error bars are Gaussian in XY ZU V W space. The resulting list of bona fide members is presented in Table 2. A list of new bona fide members confirmed in this work is given in Table 6, and their positions on the sky are displayed in Figure 3.7.

KINEMATIC MODELS OF YOUNG ASSOCIATIONS
The bona fide members compiled in Table 2 were used to build the kinematic models of young associations considered in BANYAN Σ. Objects with total error bars on their Galactic position ( σ 2 X + σ 2 Y + σ 2 Z ) above 20% of their distance or on their space velocity ( σ 2 U + σ 2 V + σ 2 W ) above 8 km s −1 were excluded. Companions in binary systems were ignored, and the spatialkinematic position of binary systems was approximated as that of the primary star to avoid the need for modeldependent mass estimates.
A rejection algorithm based on minimum spanning trees (MSTs; e.g. see Allison et al. 2009;Gagné et al. 2015b) 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 artifically 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 XY Z or U V W 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 XY ZU V W axes.
For each association containing a total of N k members, the algorithm of Cartwright & Whitworth (2004) was used to build N k +1 MSTs, in both XY Z and U V W spaces separately. The first spatial and kinematic MSTs with respective total lengths L spa and L kin include all members, and the N k additional MSTs ignore one member at a time, and have respective total lengths L i,spa and L i,kin .  Note-The average Galactic positions X , Y and Z and space velocities U , V and W correspond to the components of theτ vector defining the center of the multivariate Gaussian model. Parameters Σ00 through Σ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 detail. The covariances of the multivariate Gaussian models are given in a FITS file format in the online-only additional material.
Relative lengths L i,rel of MSTs ignoring one star each were then calculated, and an arbitrary threshold L t was set at a value 5% smaller than the 90% percentile value L i,rel 90% of all the relative MST lengths: Any individual MST with relative length L i,rel < L t 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 t = 0.9 L i,rel 90% 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 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 5.
The bona fide members that were rejected from the model construction are listed in Table 10. An average of 3 objects (typically less than 8 objects) were rejected in each young association. Only the PLE and HYA had more rejected members (10 and 15, respetively), but given their large number of members this represents less than 8% of their 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 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 ,avg and covariances C ij of all spatial-kinematic coordinates {Q i } were calculated, where the weight of a given star is set to the inverse square of its spatial and kinematic error bars σ k,spa and σ k,kin relative to the association averages σ avg,spa and σ avg,kin , added in quadrature: 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 XY ZU V W measurements bearing too much weight in the kinematic models.
The covariance matrixΣ and center vectorτ of an association were then built from its components Q i,avg and C ij , and the covariance matrix was regularized 10 to avoid numerical problems. This was done through a singular value decomposition of the covariance matrix: whereΣ sv is a diagonal matrix containing the singular values. If the determinants Σ U or Σ V were found to be negative, random noise with a standard deviation equal to half the error bars was added to the XY ZU V W coordinates of all members untilΣ was found to be nonsingular. If an association contains less than 30 members, the three spatial and three kinematic singular values are forced to a minimum of 1 pc and 0.2 km s −1 , 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 Σ ij and Σ ji to the average (Σ ij + Σ ji ) /2. In a final step, the non-diagonal elements Σ ij were forced to values within the range ± Σ ii Σ 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 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 7. The off-diagonal elements of the covariance matrices are provided in a FITS file with the BANYAN Σ algorithm as online-only material. 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: where σ 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 greenshaded 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). Modelling 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 5 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 3(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 3(d), the characteristic spatial and kinematic scales S spa and S kin are displayed for different associations as a function of age, with: where Σ spa and Σ kin are 3×3 matrices that contain only the purely spatial or kinematic terms of the covariance matrix, respectively. The values of S spa and S 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 (βPMG, 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.

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(Robin et al. , 2003(Robin et al. , 2012(Robin et al. , 2014(Robin et al. , 2017 11 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 11 Available at http://model2016.obs-besancon.fr/ thin disk population, which is based on their Model B (see their Table 5). In summary, thin disk stars are generated from a 3-slopes initial mass function and a decreasing star formation rate, and follow the evolutionary tracks of Bertelli et al. (1994Bertelli et al. ( , 2008Bertelli et al. ( , 2009 for masses larger than 0.7M , 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 semi-major axis distributions, as described by Arenou (2011) 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 modelling the field hypothesis with a mixture of N multivariate Gaussian distributions: which yields the solution described in Equation (12) for the probability P({O i }|H j,field ) associated with field component j. The resulting field probability will then be: The second problem of the approximately uniform X and Y distributions can be mitigated by artificially inflating the two diagonal elements of all covariance ma-tricesΣ 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 XY ZU V W coordinates of field stars compared to their variances (a) Distance versus age (b) Radial velocity versus age (c) Spatial size versus age (d) Velocity dispersion versus age 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 detail.
(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Σ j can be assumed diagonal and the fitting can be done simultaneously in four onedimensional spaces Z, U , V and W instead of a single 4-dimensional space. The X and Y 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 leastsquares 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 χ 2 as a function of the number of mixture components N is displayed for each dimension, and for the global fit across ZU V W . This figure shows that the goodness-of-fit does not improve significantly at N > 7 for the kinematic dimensions U V W , 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 Σ 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.
The X and Y components of the field model were arbitrarily set to a large characteristic width of 1 500 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 9 synthetic stars to be drawn and yielded a ratio of 1.69 × 10 −5 , which was divided to the total number of objects within 300 pc in the Besançon model (7.15 × 10 6 ) 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é et al. 2018b,a).

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 5 thresholds on the Galactic latitude were selected uniformly in the ranges 5-300 mas yr −1 and 0-40°, respectively. For each combination of thresholds, the ZU V W coordinates of 10 7 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 −1 ) 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 −1 ) 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. 2015b) 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.

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 as- sociation 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 k 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 simpler. 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 XY ZU V W coordinates of 10 7 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 −1 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 groups H j , j = k: where all initial priors P (H k ) and P (H field ) were set to unity. A set of 10 5 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 }) ≥ 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: where x is an array of 10 5 uniformly distributed values in the range [−1, 1]. This produces an array of thresholds where P t ≈ 0 and P t ≈ 1 are especially well sampled as ω takes larger values. Recovery rates R k (also called 'true positive rates' or TPRs) are determined as a function of threshold P t , with: where N TP is the number of synthetic stars originating from the model of association association H k with P k ≥ 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 α k that ensures P t,crit = 90% when included to the Bayesian prior is then determined: To avoid biasing the young association probabilities of ambiguous candidate members, an average of each young association factor α k , weighted by the individual membership probabilities P (H k |{O i }), is used to determine the field prior: where the sum k excludes k = field. The quantities α k are fixed for all candidate stars, but P (H 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 α k values are listed in Table 1.

THE PERFORMANCE OF BANYAN Σ AS A BAYESIAN CLASSIFIER
This section describes the classification performance of BANYAN Σ (in normal mode and in U V W -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 , 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 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.   Figure 6. Unidimensional projections of the best-fitting ten-components multivariate Gaussians 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 7 Gaussian components are displayed in blue dashed lines. Space velocity distributions are slightly asymmetric and the distribution in 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 detail.
The 7.15 × 10 6 objects within 300 pc from the Besançon Galactic model and the bona fide members listed in Table 2 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: or P k ≥ 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 the particular thresh- Figure 7. Reduced χ 2 of the multivariate field model as a function of the number of included multivariate Gaussian components. The distribution in space velocities U V W are slightly asymmetric and require 7 Gaussian components to be properly modelled, 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 detail.
olds 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 Σ 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 Mvalue is the main reason they do not perform as well as the BANYAN tools under this metric. BANYAN Σ in U V W -only mode 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 12(a) for each configuration of BANYAN Σ. In Figure 12(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 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 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,c) and the search for red L dwarfs of Schneider et al. (2017) are displayed with a 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 detail.
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 (see Table 1).
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 = 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 U V W -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 σ metrics taken in sky position and proper motion space to probabilities directly. In particular, most of the LACEwING cross-contamination is due to   . Receiver operating characteristic curve for BANYAN Σ 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 detail. Figure 11. True positive rate (TPR) as a function of Bayesian probability threshold Pt for different young associations, when only sky position and proper motion are considered. Bayesian priors were chosen for Pt = 90% to yield TPR = 50% in this scenario (e.g., see thick gray lines). The dashed thick gray lines represent the range in Pt for which the range of recovery rates are reported in Table 9. See Section 8 for more detail.
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 10. 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 constraint both XY Z and U V W sets of coordinates. Radial velocity only helps to constrain U V W . This observation was already made for BANYAN II and LACEwING (Gagné et al. 2014;Riedel et al. 2017b). 15025928-3238357 UCL 96.9 Table 8 continued Note-See Section 9 for more detail.
The models of young associations were also used to draw 10 5 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 11. This figure illustrates the effect of the α 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 (MCC; 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  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 Σ. 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 U V W kinematics that differ significantly from those of most field stars and other young associations (see Table 1). See Section 8 for more detail.
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 5 , 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 6 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 9, as the range of possible values within P t ∈ 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 9.
We investigated the effect of ignoring covariances between the input observables (see Section 3.3) with a 10 4elements Monte Carlo simulation. The median proper motion and parallax errors of all bona fide members in Gaia-DR1 (respectively 0.2 and 0.15 mas yr −1 , 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 were attributed 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 4 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 −5 % in probability, 3 × 10 −6 % in optimal distance and 8 × 10 −6 % 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.

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 were 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. (2017a) 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 5 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 8 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. (2017a), 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.

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. It is also made publicly available in IDL and Python at www.exoplanetes.umontreal.ca/banyan/ banyansigma.php. Additional figures and information on this work can be found on the website www.astro. umontreal.ca/~gagne.
We strongly encourage users to investigate the position of candidate members in appropriate colormagnitude 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 10). 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é, J. K. Faherty,É. Artigau et al., in preparation; see also Gagné et al. 2017b andGagné et al. 2015a for preliminary results from BASS-Ultracool), as well as a search for new stellar members based on Gaia-DR1 (J. Gagné, O. Loubier, J. K. Faherty et al., 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. Note-True-positive rates (TPRs), false-positive rates (FPRs) and Matthews correlation coefficients (MCCs) are reported at the 'critical' Pt =90% threshold, and their reported ranges are for thresholds in the range 85-95%. See Section 8 for more detail.
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. (2015c), 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;Kirkpatrick et al. 2003), 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). EEM acknowledges support from the NASA NExSS program.
JG designed BANYAN Σ, compiled the bona fide members, wrote the IDL codes, wrote the manuscript, generated figures and led all analyses;EEM shared parts of the young association literature data, characteristics and bona fide members, and provided general comments;ORL wrote the initial Python translation of BANYAN Σ; LM, AR and DR performed the BANYAN I, LACEwING and convergent point tool calculations used in Section 8, respectively; AR also provided help with the bona fide members compilation; DL and JKF shared ideas and comments; LP shared ideas and provided comments especially for the ROC curves analysis; AR performed custom Besançon Galactic simulations and wrote the second pagraph of Section 6; and RD shared comments and supervized ORL.
Error bars on the optimal radial velocity σ ν and distance σ can be defined by measuring the characteristic width of P o ({O i }|H) along ν and in the vicinity of (ν o , o ). The effect of the Jacobian term 4 will be ignored to obtain an analytical approximation of (σ ν , σ ).
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 σ ν and σ : In the case of σ ν , this yields: leading to: The case of σ requires the introduction of a new variable that is defined in the range ] − ∞, ∞[ and matches the distance = for ≥ 0. Assuming that o >> σ will ensure that the Bayesian likelihood P 0 (ν, ) ≈ 0 for all negative values of . It follows that: leading to: 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.