# Untangling the Galaxy III: Photometric Search for Pre-main Sequence Stars
with Deep Learning

Aidan McBride Department of Physics and Astronomy, Western Washington
University, 516 High St, Bellingham, WA 98225  Ryan Lingg Department of
Computer Science, Western Washington University, 516 High St, Bellingham, WA
98225  Marina Kounkel Department of Physics and Astronomy, Western Washington
University, 516 High St, Bellingham, WA 98225  Kevin Covey Department of
Physics and Astronomy, Western Washington University, 516 High St, Bellingham,
WA 98225  Brian Hutchinson Department of Computer Science, Western Washington
University, 516 High St, Bellingham, WA 98225  Computing and Analytics
Division, Pacific Northwest National Laboratory, 902 Battelle Blvd, Richland,
WA 99354<EMAIL_ADDRESS><EMAIL_ADDRESS>
<EMAIL_ADDRESS>
###### Abstract

A reliable census of of pre-main sequence stars with known ages is critical to
our understanding of early stellar evolution, but historically there has been
difficulty in separating such stars from the field. We present a trained
neural network model, Sagitta, that relies on Gaia DR2 and 2MASS photometry to
identify pre-main sequence stars and to derive their age estimates. Our model
successfully recovers populations and stellar properties associated with known
star forming regions up to five kpc. Furthermore, it allows for a detailed
look at the star-forming history of the solar neighborhood, particularly at
age ranges to which we were not previously sensitive. In particular, we
observe several bubbles in the distribution of stars, the most notable of
which is a ring of stars associated with the Local Bubble, which may have
common origins with the Gould’s Belt.

††software: TOPCAT (Taylor, 2005), Pytorch (Paszke et al., 2017)

##  1 Introduction

Historically, the pre-main sequence (PMS) stars that have been the easiest to
identify and classify are those that are the youngest and are still in
possession of their natal envelopes and/or protoplanetary disks. Such sources
could be identified on the basis of large infrared excess, and these dusty
young stellar objects (YSOs) have been searched for using a number of all-sky
surveys, such as using IRAS, 2MASS, AKARI, and WISE, (e.g., Prusti et al.,
1992; Koenig et al., 2012; Tóth et al., 2014; Marton et al., 2016).
Furthermore, detailed infrared maps of a large number of star forming regions
have been constructed with more targeted surveys, such as with Spitzer and
Herschel (e.g., Evans et al., 2009; Megeath et al., 2012; Fischer et al.,
2017).

However, after a star loses its protoplanetary disk, its colors begin to
resemble those of much more evolved field stars, making follow-up
identification difficult.

###  1.1 Gaia DR2 classification of YSOs

In comparison to previously available techniques (such as spectroscopic follow
up measurements, or X-ray emission), the release of Gaia DR2 (Gaia
Collaboration et al., 2018) has allowed a revolution in the search and
characterization of young stars. Through its unprecedented precision in the
measurements of parallax and proper motions, as well as its remarkable
photometric quality, two new techniques became available to the community.
First is the phase space clustering. Young stars form in the dynamically cold
molecular clouds. These clouds would typically form anywhere from a few
hundred to several thousands of stars in relatively close proximity that would
typically have low velocity dispersion. Thus, through searching for an
overdensity in the position and velocity space it is possible to identify a
young comoving group of stars. Such clustering has been employed both
systematically across the entire Galactic Disk (Kounkel & Covey, 2019,
hereafter Paper I) as well as to better constrain the membership of individual
star forming regions (e.g., Kounkel et al., 2018; Galli et al., 2019; Damiani
et al., 2019).

But, clustering requires that all of the stars in a comoving group retain the
group’s characteristic velocity in order to be identifiable. As these groups
slowly dissolve into the Galaxy and lose coherence, an increasingly small
fraction of them is recoverable. Indeed, even 1 Myr populations have some
stars that have already been ejected from the clusters they inhabit (McBride &
Kounkel, 2019; Schoettler et al., 2020; Farias et al., 2020). Searching for
such high velocity YSOs may be of interest to better characterize intracluster
dynamics, but it is impossible to do through clustering. Furthermore, some
young populations may be too small or too diffuse to robustly identify them
with clustering at all.

The second method that Gaia DR2 made possible is through examining the
position of the stars on the HR diagram. YSOs are overluminous compared to
their main sequence counterparts due to their still-inflated radii, and most
are fainter and cooler than the red giants. If the distance is known
accurately, it is possible to resolve the degeneracy between the nearby dwarfs
and distant giants, and thus separate YSOs from more evolved stars.

Such an approach is rather simple to use when attempting to identify YSOs in a
single star forming region with a known position on the sky and known
distance, particularly if this region is only a few Myr old. In this case, the
low mass PMS stars can be cleanly separated from the low mass main sequence
counterparts, and it is possible to determine color cuts that would prevent
contamination from red giant stars and or high mass main sequence stars.
However, extending this to multiple populations that have different ages,
distances, or extinctions is difficult. Zari et al. (2018) performed a
selection of YSOs in the Gaia DR2 catalog that are consistent with being
younger than the 20 Myr isochrone and that are located within 500 pc. Although
the catalog effectively identifies sources throughout known nearby star
forming regions, at larger distances the contamination does become
significant. Thus, it is necessary to reevaluate the selection criteria if one
wishes to reliably extend the catalog beyond 500 pc.

Machine learning, and, in particular, the use of neural networks, is a method
that facilitates the search for complex correlations in large volumes of data.
Machine learning classifiers have been used to search for young stars in a
number of works, from searching for infrared excess (Marton et al., 2016,
2019; Chiu et al., 2020), to using H$\alpha$ in conjunction with photometry to
search for Herbig Ae/Be stars (Vioque et al., 2020), to using the optical
Hubble Space Telescope colors to give a probabilistic assessment of young
stars in the Magellanic clouds (Ksoll et al., 2018).

###  1.2 Derivation of stellar ages

Beyond classifying a star as young, extracting its properties (such as its
age) can be a challenge. The way this is commonly done is through comparison
of photometric colors (or age-sensitive spectroscopic features) to theoretical
isochrones. While this practice has a long standing history (e.g., Cohen &
Kuhi, 1979; Greene & Meyer, 1995; Covey et al., 2010; Da Rio et al., 2012),
this process is not trivial (e.g., Olney et al., 2020), especially with the
inconsistencies between young stars discussed earlier.

The first issue lies with the isochrones themselves. Over the years, a number
of different stellar evolution models have been developed (e.g., D’Antona &
Mazzitelli, 1994; Baraffe et al., 1998; Siess et al., 2000; Baraffe et al.,
2015; Choi et al., 2016), and the ones that seem to have gained the most wide-
spread usage in the community in the recent years are the PARSEC isorchrones
(Marigo et al., 2017). Due to slightly different assumptions regarding the
underlying stellar physics, these models produce distinct isochrones and
evolutionary tracks, and thus produce different age and mass estimates even
when applied to the same stellar population (Hillenbrand et al., 2008).
However, no isochrones offer a perfect match to the data, especially for the
low mass stars.

For example, the M dwarfs appear to be overinflated compared to what the
isochrones would suggest even in one of the best studied open clusters, the
Pleiades (Jackson et al., 2018). Thus, attempting to estimate their age
through isochrone fitting would yield a systematically younger age than what
is appropriate for the cluster. Similarly, a cluster “birthline” (i.e., the
region of the parameter space that would correspond to a 0 Myr population) is
ill-defined, such that in the young populations that are just a few Myr old,
the higher mass stars appear to be systematically older than their low mass
counterparts (Hartmann, 2009; Herczeg & Hillenbrand, 2015). The presence of a
protoplanetary disk further alters the photometry in such a way that makes it
difficult to place a star onto the isochrones correctly.

The second issue lies with the physical properties of YSOs. They tend to be
more complex than even the most advanced stellar evolution models can account
for. Due to being mostly convective, a large fraction of their photospheres
are covered with spots, resulting in a mismatch in the effective temperature
$T_{eff}$ and the expected mass of the star. Finally (albeit unavoidably) is
an issue of stellar multiplicity. Among the main sequence stars, the binary
sequence can be clearly seen and separated from the single stars. However, in
young populations with an intrinsic age spread of even a few Myr, the visual
binaries may be inferred to be systematically younger than what is appropriate
(e.g., Bouma et al., 2020).

Third, there is an issue of the self-consistency of the fitting process. Even
for the same populations of stars, using the same set of isochrones, but
focusing on somewhat different features and using a different interpolation
method, it is possible to produce an age estimate that is somewhat
inconsistent between various works. Taking Orion as an example, particularly
the region near 25 Ori cluster, while different authors were able to estimate
roughly comparable ages (Kounkel et al., 2018; Briceño et al., 2019; Zari et
al., 2019), some infer ages that are systematically older (Kos et al., 2019).
Such differences would only compound when comparing ages of individual stars.

To compensate for some of these issues, data driven models may perform better
compared to the theoretical isochrone fitting. With distilling the previously
existing estimates of ages for stars (both in the cases when the age can be
assigned to all stars in a cluster on a population level, and in cases where
measurements for individual stars are available), it is possible to construct
a neural network that would assign ages to stars. While the predictions it
would generate could only be as accurate as input data on which the network is
trained on, through leveraging the ages derived by various methods, the
systematic differences can be significantly reduced. This includes systematic
differences between low and high mass stars, as well as differences between
various stellar evolution models, resulting in a more self-consistent
interpolation that is more faithful to the data.

In this work, we present a tiered deep learning model, that we refer to as
Sagitta. This model identifies the PMS stars using Gaia DR2 and 2 Micron All-
Sky Survey (2MASS) photometry and astrometry and estimates the ages of these
stars. In Section 2 we describe the data that were used to train Sagitta, as
well as the data on which we perform the evaluation. In Section 3 we detail
the process of constructing and training the model. In Section 4 we test the
results benchmarked against other catalogs of PMS stars, as well as known star
forming regions. In Section 5 we discuss the features in the data, such as
their implications on the star forming history of the solar neighborhood and
the origin of the Gould’s belt. Finally, we conclude in Section 6.

##  2 Data

Neural networks are reliant on labelled training data in order to generate a
model able to generalize; that is, make accurate predictions for new sources.
The bulk of the training data was obtained from Kounkel et al. (2020,
hereafter Paper II). This catalog contains almost 1 million stars that have
been clustered together into more than 8000 different moving groups, extending
up to parallax limit of $\pi>0.2$ mas. Average ages (ranging from $<1$ Myr to
$1$ Gyr) has been inferred through the photometry of each moving group’s
members using the Auriga neural network, which is described in full in Paper
II.

To train its sister network, Sagitta, that we present in this paper, we ingest
8 parameters for each star in the sample. Six of them are the photometries in
different bandpasses, namely $G$, $G_{BP}$, $G_{RP}$, $J$, $H$, and $K_{s}$,
with the 2MASS passbands queried using the precomputed Gaia DR2 crossmatch
table. We also rely on parallax, as well as an approximation of extinction
$A_{V}$. The latter is not necessarily the intrinsic extinction the star may
have (e.g., it does not account for the presence of a protoplanetary disk),
but rather, it is an estimate along the line of sight to the star based on its
galactic positions $l$, $b$, and $\pi$ (Section 3.3).

The tasks of both classification (i.e., blind selection of PMS stars) and
regression (i.e., interpolation of the ages for the identified young stars)
are reliant on the different features of this catalog; furthermore, they
require different augmentation procedures to improve the homogeneity of
coverage.

###  2.1 Classification training sample

The first necessary step is to select the PMS stars in the catalog from Paper
II. Massive YSOs will reach the main sequence sooner than the low mass ones:
e.g., OB stars will be born directly on the main sequence, whereas M dwarfs
may take as long as 100 Myr to reach it. As main sequence stars of similar
mass are difficult to distinguish from one another, regardless of their age, a
simple cut in age is not sufficient to reliably separate PMS stars. Rather, we
compared the HR diagrams of populations in each 0.1 dex age bin to that of the
Pleiades to determine the most massive/bluest star that can still be
considered pre-main sequence. We then assign the extinction-corrected
$[G_{BP}-G_{RP}]_{0}$ color corresponding to such star as the location of the
turn-off for that age bin. The extinction was estimated on a population level
in Paper II. We then interpolated across all the age bins to obtain the cut of

$[G_{BP}-G_{RP}]_{0}>49.3686-14.3347\times t+1.05042\times t^{2}$

where $t$ is the age of the population in dex. This relation is valid for
$t<7.85$ dex (Figure 1). Using this age-dependent relation to identify the
color at which the MS-PMS transition should occurs at a given age, we assign a
preliminary PMS designation to all low-mass stars redward of the critical
color in each population, according to its age.

Figure 1: A demonstration of a color-dependent cut-off for pre-main sequence
stars as a function of age in comparison to the Pleiades. Note the presence of
the binary sequence in each of the populations.

Furthermore, we imposed a cut of

$M_{G_{0}}>2.8\times[G_{BP}-G_{RP}]_{0}-2$

to separate possible contamination from red giants.

A total of 62,484 sources (about 6%) of the catalog in Paper II meet these
criteria, and are identified as the training sample of likely PMS stars for
Sagitta. The remaining sources are considered to be evolved, i.e., main
sequence or post-main sequence sources. All sources receive a binary numerical
flag to indicate their evolutionary state, with 1 assigned to all sources
satisfying the YSO criteria and 0 assigned to the remaining evolved sources.

While the catalog from Paper II offers a comprehensive coverage of PMS stars,
it is incomplete at the oldest age bins for stars that are more evolved. All
moving groups eventually dissolve into the Galaxy, and after this happens,
they are no longer recoverable through clustering. Therefore, the distribution
and properties of the red giants that are found in the field is not
represented by the stars in the Paper II catalog. This introduces a bias that
may result in reddened red giants that scatter into PMS parameter space (or
even located well above it) to erroneously be classified as PMS, as the
classifier would not know enough about red giants to reject them. To correct
for this bias we took 3 million randomly selected stars from Gaia DR2 (using
the random_index), with the same quality cuts as the ones specified in Section
2.3 to better train the network on the colors of older field stars to be able
to discriminate against them. Any stars that happened to also be included in
the clustered catalog were excluded from this random selection. We assume that
the remaining stars in this random catalog can all be considered to be
evolved. If any true PMS stars remained in this random sample, their fraction,
as well as the overall number is expected to be so small as not to make a
substantial difference for the classifier, as PMS are comparatively rare
relative to the older stellar classes.

The spatial position of the stars (i.e., $l$ and $b$) was not used directly in
the training sample to avoid introducing spatial bias. However, some positions
can be inferred through the combination of $\pi$ and $A_{V}$. To prevent
spatial bias, we used augmentation, i.e., the process of making the training
catalog larger through artificially modifying it to reduce various biases.
Augmentation is also beneficial to improve performance beyond 1 kpc, where the
training sample is highly incomplete, as few PMS stars at those distances are
found in the Paper II catalog due to the sensitivity limits. To augment the
sample, we modified the photometry of the real stars (both the clustered and
random catalogs) to simulate the effects of them being found at different
distances (up to 5 kpc) and with different $A_{V}$ (up to 10 mag), both of
which were randomly generated, and random errors in flux drawn from a normal
distribution were added to all passbands. The extinction coefficients for
different passbands were taken from the web interface of the PARSEC isochrones
(Marigo et al., 2017) as they cover a wide range of stellar masses and ages,
with various passbands applied to the photometry. Each real star was drawn
multiple times, with this multiplier treated as a hyperparameter in the model
for each subset (Section 3.4). These ‘synthetic’ stars were passed along to
the classifier alongside with the photometry of the real stars to improve
generalization.

###  2.2 Regression training sample

The training sample for the regression to estimate ages was limited only to
the stars identified as PMS stars in the catalog from Paper II from above,
excluding the sources considered to be evolved. Otherwise, since the main
sequence stars of similar mass share the parameter space in fluxes, retaining
them in the sample could introduce considerable biases to the PMS stars.

Most of the PMS stars in Paper I and Paper II are found in stellar strings,
i.e., extended populations spanning several tens or even hundreds of pc. Each
individual region can sustain star formation for up to $\sim$10 Myr. Such a
duration is hardly noticeable within the errors of the ages assigned to
populations older than 100 Myr, i.e, the colors and fluxes of a 90 Myr and 100
Myr stars are not going to be very different. However, the differences in
fluxes between the youngest and the oldest generation of stars in the same
region is more pronounced in regions that are still forming, that is to say, a
1 Myr and 10 Myr old stars look very different. Therefore, assigning a single
age to all the stars in a single string if there is any sort of underlying age
gradient (as is the case in Orion or Sco Cen, for example) can introduce some
biases. In estimating ages for the moving groups in Paper II, Auriga
preferentially considers the oldest stars in a region, overestimating the ages
of the younger stars. Without correcting for this effect in the training
sample, Sagitta cannot accurately estimate ages of stars younger than a few
Myr.

Thus, to compensate, several steps were taken. First, many of the strings can
be subdivided into populations more homogeneous in age. In Paper I, some of
the strings had to be manually assembled from smaller subgroups based on the
coherence in phase space. Furthermore, the HR diagrams of them were visually
examined in trying to identify populations of different age sequences in close
proximity (e.g., as is the case with $\rho$ Oph and Upper Sco), and some
attempt was made to separate them into subgroups. We used Auriga on these
subgroups to generate a somewhat more granular distribution of ages in the
training sample.

Table 1: Ages different from Paper II assigned to young populations in the training set.  Region | Source | Age (Myr)
---|---|---
Ara OB1a | Wolk et al. (2008) | 3
Carina: Tr 16 | Smith & Brooks (2008) | 3
Cep OB2a | Kun et al. (2008) | 7
Cep OB2b | Kun et al. (2008) | 3.7
Cep OB3b | Kun et al. (2008) | 4
Cep OB6 | Kun et al. (2008) | 38
Chamaeleon | Luhman (2008) | 2
CrA | Neuhäuser & Forbrich (2008) | 6
Cyg OB1 | Reipurth & Schneider (2008) | 7.5
Cyg OB2 | Reipurth & Schneider (2008) | 5
Cyg OB3 | Reipurth & Schneider (2008) | 8.3
IC 348 | Bally (2008) | 2
IC 1396 | Walawender et al. (2008) | 1
IC 5146 & W4 | Herbig & Reipurth (2008) | 1
LK H$\alpha$ 101 | Andrews & Wolk (2008) | 0.5
Lagoon Nebula | Tothill et al. (2008) | 1
Lower Cen/Crux | Preibisch & Mamajek (2008) | 16
Lupus | Comerón (2008) | 3.2
Monoceros | Carpenter & Hodapp (2008) | 6
NGC 1333 | Walawender et al. (2008) | 1
NGC 2264 | Dahm (2008) | 3
NGC 6383 | Rauw & De Becker (2008) | 2
NGC 6604 | Reipurth (2008a) | 4.5
NGC 6823 | Prato et al. (2008) | 5
Per OB2 | Bally (2008) | 6
Rosette Nebula | Román-Zúñiga & Lada (2008) | 3
Serpens | Herczeg et al. (2019) | 3
Sh 2-234-Stock 8 | Reipurth & Yan (2008) | 2
Taurus/Auriga | Kenyon et al. (2008) | 1
Upper Cen/Lup | Preibisch & Mamajek (2008) | 17
Upper Sco | Preibisch & Mamajek (2008) | 5
$\rho$ Oph | Wilking et al. (2008) | 0.3

Second, to achieve a better consistency with the ages of well-studied star
forming regions, we identified the moving groups from Paper II that correspond
to the populations listed in the Handbook of Star Forming Regions (Reipurth,
2008b, c), and assigned them the more appropriate ages that are reported in
the Handbook (see Table 1 for group identifiers and corresponding ages).

Third, individual ages of some young stars are available in the literature.
Namely, we’ve included the sources in the catalogs from Palla & Stahler
(2002); Kun et al. (2009); Delgado et al. (2011); López Martí et al. (2013);
Fang et al. (2013); Kumar et al. (2014); Herczeg & Hillenbrand (2014); Getman
et al. (2014); Erickson et al. (2015); Azimlu et al. (2015); Fang et al.
(2017); Suárez et al. (2017); Prisinzano et al. (2018); Panwar et al. (2018)
that have reliable parallaxes and that meet the age-dependent criteria for a
source to be identified as PMS from Section 2.1. This added 6,248 stars.

Finally, the ages have been reevaluated in the Orion Complex. As this region
singularly contains the largest number of stars out of any other population in
the catalog, and it contains stars that span in age from $<$1 to 12 Myr, it is
of particular importance for the training sample from Orion to achieve good
accuracy across different age bins. The region does have a rather complicated
morphology that could not be broken into subgroups in Paper I. However, a
different, more granular analysis was performed by Kounkel et al. (2018),
using hierarchical clustering of the 6-dimensional phase space to segment the
Complex into 190 different groups, and an average age was estimated for each
group. The sample in that work, however, is limited to the stars almost a
magnitude brighter than the sample in Paper II. Excluding the fainter stars
does introduce a bias in the training process. Thus, to shuffle these low mass
stars into the most appropriate group, we created a simple fully connected
neural network that has one layers with 300 neurons, taking in $\alpha$,
$\delta$, $\mu_{\alpha}$, $\mu_{\delta}$, and $\pi$, and outputting a
probability of belonging to each one of the 190 groups. The members of the
Orion Complex from Paper II that were not included in the catalog from Kounkel
et al. (2018) were then assigned to the group with the highest probability.
Then, each star was given a label of the average age of the group it was in.

Similarly to the classification training sample, to help negating some of the
potential spatial biases, the training sample for the age regression was
augmented. As above, the synthetic sample was generated from the real data
through simulating the effects of changing the distance and extinction.

###  2.3 Evaluation sample

To evaluate the performance of the classification and regression models, we
downloaded the Gaia DR2 data that satisfied the following quality criteria:

$\texttt{phot\\_bp\\_rp\\_excess\\_factor}>1+0.015\times\texttt{bp\\_rp}^{2}$
$\texttt{phot\\_bp\\_rp\\_excess\\_factor}<1.3+0.06\times\texttt{bp\\_rp}^{2}$
$\texttt{ruwe}<1.4$ $\texttt{phot\\_g\\_mean\\_flux\\_over\\_error}>10$
$\texttt{phot\\_bp\\_mean\\_flux\\_over\\_error}>10$
$\texttt{phot\\_rp\\_mean\\_flux\\_over\\_error}>10$ $\texttt{parallax}>0.2$
$\texttt{parallax}/\texttt{parallax\\_error}>10~{}\textrm{or}~{}\texttt{parallax\\_error}<0.1$

These criteria are adapted from the quality cuts by Lindegren et al. (2018).
They also ensure that the sources are nearby enough to have a complete
coverage of the volume of space over which low mass PMS stars are detectable.
All three fluxes from Gaia DR2 are required to be detected with high signal to
noise. 2MASS fluxes can be undetected, in which case they are set to a
constant outside of their maximum range (See Section 3.2).

The resulting sample consists of $\sim$139.3 million stars. We note that the
models presented in the paper can be expected to work even if some of the
selection restrictions are relaxed, although it is not done in this paper for
the purposes of the computational expediency.

Figure 2: Top: The distribution of PMS stars on the sky that were used as part
of the training process. The sources are color coded by the age assigned to
them. Bottom: HR diagram of stars in the input dataset, uncorrected (left) and
corrected (right) for extinction. Grey points represent evolved stars used in
the classifier training, while color coded points represent pre-main sequence
stars weighted by isochrone age, used to train the age regressor model.

##  3 Sagitta neural network

###  3.1 Network Architecture

Figure 3: Architecture of the model used by Sagitta. This architecture is used
by each of three networks in this work, independently of each other.

For flexibility and effectiveness in model structure, we implement a set of
three convolutional neural networks (CNNs) where each network serves a
distinct purpose. The first network generates the extinction map based on
galactic position. The second network assigns each star a probability of being
pre-main sequence based on its photometry (using the training sample described
in Section 2.1). And the third network approximates the ages of PMS stars,
typically ranging from $<1$ Myr to $>40$ Myr (using the training sample
described in Section 2.2).

All three networks within Sagitta share a common architecture adapted from the
MNIST configuration of an online repository of
models111https://github.com/eladhoffer/convNet.pytorch/blob/master/models/mnist.py,
which was further adapted in Paper II for Auriga. However, there are important
differences between these models. The aforementioned examples use
2-dimensional arrays as inputs (e.g., for Auriga the inputs are the the fluxes
across various bands for each star in a cluster). In the case of Sagitta
though, the input is only a 1-dimensional array representing parameters of a
single star, so the operational dimensions (such as convolution and pooling)
in the network have been reduced.

The networks were implemented using PyTorch (Paszke et al., 2017). The model’s
architecture broadly consists of three main segments (Figure 3). The first
segment is made up of of four 1-dimensional convolutional layers that
gradually increase the number of channels to a maximum size of 128. In this
segment, only the first and third convolutional layers are followed by max-
pooling layers, but each of the convolutional layers or convolutional and max-
pooling layer pairs is followed by ReLU and Batch-Normalization operations.
The second segment consists of a convolution layer that reduces the number of
channels down from the output of the first segment followed by a data
reshaping operation that transforms the convolution layer’s output into a
1-dimensional list (width $\times$ channels $\rightarrow$ width only). Once
reshaped, the data are then fed into the third segment which consists of three
fully connected linear layers where each layer is followed by a ReLU
operation. These final layers expand and then contract the processed data down
to a single scalar output value. For the $A_{V}$ estimation and age estimation
models the output value is kept as is (i.e., linear output activation),
however in the case of the YSO classifier model this final value is then fed
into a logistical sigmoid function to bound the output probabilities between 0
and 1.

###  3.2 Data Handling

In the process of passing the data through the model, it is beneficial to
first normalize all of the input parameters to a similar range with mean close
to zero. This creates a more comparable dispersion of input parameter
magnitudes and mitigates potential issues with numerical stability or
inherently biasing any input because of its original scaling. Thus, all the
parameters for both classification and regression were linearly scaled to the
range of $[-1,1]$ based on the lower and upper bounds specified in the Table
2. Although the overall distribution of either parallaxes or fluxes is not
Gaussian, thus resulting in a skewed distribution that is dominated by values
towards one of the ends of the normalization, the overall bounds were
considered to be sufficiently effective.

Table 2: Normalization constants used in training.  Parameter | Lower | Upper
---|---|---
$\pi$ (mas) | 0 | 29
$A_{V}$ (mag) | 0 | 5
$G$ (mag) | 4 | 20
$G_{BP}$ (mag) | 0 | 21
$G_{RP}$ (mag) | 0 | 19
$J$ (mag) | 0 | 18
$H$ (mag) | 0 | 18
$K$ (mag) | 0 | 18
Age (dex) | 6 | 8

A number of sources in the catalog may have one or more fluxes missing (most
commonly in the 2MASS bands). The non-detection usually carries meaningful
information, as it shows that the source is fainter than the detection limit,
and limiting the catalog to only the sources for which complete data across
the different bands are available would not be optimal. Nonetheless, neural
networks are unable to handle null values as inputs.

Thus, to allow the inclusion of sources with incomplete data, we set the
missing values of both the training and evaluation set to the upper limit
specified in Table 2. These upper limits are typically somewhat fainter than
the detection limit in each band, allowing the network to learn to give these
fluxes an appropriate weight.

Another aspect that was considered when structuring the data was the exact
ordering of the stellar input parameters given to the YSO Classifier and Age
Regressor networks. Convolutional layers in a CNN operate with the use of a
sliding filter component, wherein only data spatially close enough together to
fit inside the filter can be used to directly detect patterns in that layer.
The input features were ordered as $A_{V}$, $\pi$, $G$, $G_{BP}$, $G_{RP}$,
$J$, $H$, and $K$ to preserve the rough ordering of the bandpasses with the
wavelength, and to keep $\pi$ and $G$ adjacent as they are measured in the
same dataset, and attaching it at the end would have resulted with $\pi$ being
associated with occasionally incomplete 2MASS data. Similarly, $A_{V}$ has the
strongest effect in the optical portion of the spectrum. However, as the
initial width of the sliding filter is 5 elements being convolved together,
the order should not have a very strong effect on the final results. Given
that the features are not spatial in the traditional sense, we also tried
using a fully connected deep neural network instead of the CNN, but found it
to perform somewhat worse.

###  3.3 $A_{V}$ estimate

In order to help with differentiating the PMS stars from reddened massive main
sequence stars and red giants that scatter into the parameter space PMS stars
inhabit, one parameter that can help is an estimate of $A_{V}$ along the line
of sight.

This estimate was obtained from the neural network used to generate the
completeness map in Paper II. The model was not modified in any way; rather,
it is ported into Sagitta directly as the first step.

The $A_{V}$ estimator uses the same architecture as the classifier and the age
regressor. It was trained on 3 million randomly chosen stars from Gaia DR2
that had the same quality constraints as the ones imposed in 2.3 and with
measured $A_{G}$ reported in the catalog.

The network used $l$, $b$, and $\pi$ to predict $A_{V}$ (scaled from $A_{G}$
by a factor of 0.859, Marigo et al., 2017) corresponding to the particular 3d
spatial position. Although transformation extinction from one bandpass to
another can be a complex process, the linear transformation was done for the
sake of nomenclature. As all of the parameters are also normalized through
linear scaling, the net result is comparable to training on $A_{G}$ directly.

In training, positions were normalized from 0 to 1 for $l$ from 0 to 360∘, $b$
from -90 to 90∘, $\pi$ from 0 to 5 mas, and $A_{V}$ from 0 to 5 mag, and the
individual measurements of $\pi$ or $A_{V}$ were allowed to exceed the maximum
to be $>1$ after the normalization. The training was done using the Adam
optimizer, mean square error loss, and a learning rate of $10^{-3}$.

The resulting estimate of $A_{V}$ is consistent to within 0.3 mag with the
extinction map from Green et al. (2019) over the applicable volume, as well as
to the cluster $A_{V}$s estimated on the population level in through pseudo-
isochrone fitting with Auriga (Paper II). The resulting 3-dimensional
extinction map is shown in the Figure 4.

This is sufficiently precise for the purpose of this paper, as both the
classifier and the age regressor do not depend on the absolute magnitude of
$A_{V}$ (or $A_{G}$) directly. Rather, they rely on the non-linear
correlations that are present in the data that can be inferred with a help of
this parameter, and they can learn to compensate for color-dependent
systematic differences that may be present.

Figure 4: A slice of sphere across the galactic plane, showing the
3-dimensional distribution of the sources in the evaluation sample, color-
coded by the predicted $A_{V}$ values.

###  3.4 YSO Classifier

####  3.4.1 Training

The classification network was trained to perform a binary classification
task, with classes 0 (not PMS) and 1 (PMS). It maps each stellar parameters
$A_{V}$, $\pi$, $G$, $G_{BP}$, $G_{RP}$, $J$, $H$, and $K$ to a single scalar
output in $(0,1)$ representing the star’s probability being PMS, using a
logistic sigmoid output activation. The model is trained to maximize the log-
likelihood of the training data (equivalent to cross-entropy loss).

The data from Paper II was partitioned into three disjoint sets where each
served a distinct purpose in classifier development. The training set,
containing the first $80\%$ of the sources, was comprised of the stars for the
model to train off of. The development set, containing the next $10\%$ of the
sources, was used during training to evaluate model performance but was never
shown to the model as examples. The testing set, containing the last $10\%$ of
sources, was only used after tuning to confirm generalization.

With the training set in place, we augmented the catalog in order to increase
performance. The subsets that were sampled for augmentation included PMS stars
from Paper II, the non-PMS stars from Paper II, and the randomly selected 3
million star sample from Gaia DR2 (see Section 2.1. Due to the difference in
size of these subsets, all the stars in each of these subsets was sampled
certain amount of times, changing $\pi$ and $A_{V}$ to affect the flux, to
produce the augmented sample. The number of times each star was sampled was
treated as a hyperparameter in the training process 3. Multiple models were
trained on all permutations of the augmentation ratios, however adjusting the
ratios did not appear to have a strong impact on the performance of the model.

Finally, the evaluation sample described in section 2.3 (that consisted mostly
of sources for which we did not have explicit a-priori labels) was used for
final testing and compare performance of various models through examining
known star forming regions, regions of high extinction, and other features of
the solar neighborhood (Section 3.4.2).

To increase the classifier’s performance in distinguishing PMS stars, during
the data augmentation phase we oversampled PMS stars (i.e., used stratified
sampling), yielding 15% PMS stars during training (compared to $1.6\%$ in the
training set). By improving class balance, the model had to focus more on
correctly identifying and reducing contamination in the PMS star class.
However, while the initial augmentation has improved sensitivity to more
distant PMS stars compared to the unaugmented sample, continuing to grow their
number in training through augmentation would not necessarily result in a
better classifier, as it provides diminishing returns. Once the training had
completed, the output detection threshold used for extracting the PMS stars
was then selected based on methods described in section 3.4.2. Fine tuning of
the PMS to non-PMS star ratio was found to not have a significant impact on
the performance of the model.

Table 3: Classifier Hyperparameter Tuning Values  Hyperparameter | Values
---|---
$Optimizer$ | Adadelta, Adagrad, Adam, RMSProp, SGD
$Learning$ $Rate$ | 0.001, 0.01, 0.1
$Dropout$ | 0%, 10%, 30%, 50%, 70%
$Minibatch$ $Size$ | 5000, 10000, 25000, 50000
$Weight$ $Decay$ | 0, 0.00001, 0.0001, 0.001
Subset | Star Sample Rate
PMS stars from Paper II | 20, 25, 50
Non-PMS stars from Paper II | 3, 5, 10
Random stars from Gaia DR2 | 1, 2

Hyperparameter tuning and early stopping were also employed to help improve
classifier predictions. The list of possible settings for each of the
hyperparameters tuned are listed in table 3. The models in the hyperparameter
sweep were trained until the development set loss failed to beat its best loss
for 20 successive epochs, at which point that model training stopped. During
each model’s training, only the snapshot of the model with the best
development set performance was saved. Once the sweep finished, each model’s
predicted outputs on the development set was used to visually confirm that the
model was predicting desired values. The final model used in the pipeline was
chosen based off its low development set loss and qualitatively good
predictions.

Through the hyperparameter sweep, it was found that using very little or no
weight decay consistently provided models with the best development set
performance. All of the other hyperparameters tuned on seemed to not produce
any significant improvement in model performance one way or another. The
configuration of hyperparameters that produced the best model was comprised of
Adagrad for the optimizer, $0\%$ dropout, a batch size of $5000$, a learning
rate of $0.01$, and a weight decay of $10^{-5}$.

####  3.4.2 Classifier validation

The Upper Sco region, which contains the $\rho$ Oph dark cloud, is
particularly useful in evaluation of how reliably the classifier can
discriminate between bona fide YSOs and the contamination from evolved stars
that have been reddened to the PMS parameter space. It is located nearby, with
$\pi\sim 7$ mas. No other star forming regions are known to be located behind
it, nor are there likely to be any distant undiscovered populations behind it,
given its high elevation above the Galactic Plane, at $b\sim 20^{\circ}$.
Therefore, if any sources identified as PMS are located far beyond, e.g., 200
pc, they are most likely to be false positives. The precision in distance that
can be inferred from the parallax decreases the further away a star is, while
the number of field stars in each parallax bin increases. Therefore, the large
parallax of Upper Sco makes it particularly easy to separate the stars
associated with this region from the false positives, to the degree that even
other star forming regions along the Gould’s Belt do not.

Figure 5:  Figure 6: Identified PMS sources $>70$% probability towards $\rho$
Oph and Upper Sco, plotted over the Gaia DR2 map of the sky (Gaia
Collaboration et al., 2018). As the $\rho$ Oph dark cloud has high extinction,
it is clearly visible in this map. Note the highest confidence PMS sources are
tracing the known regions of star formation. On the other hand, sources with
lower probability tend to be co-located along the pertruding filaments that
are not actively forming stars.

Figure 7: HR diagram of the evaluation sample. Top: Color coded by the maximum
probability of identifying a star as PMS in each 2-dimensional bin, in two
distance slices. Bottom: Color coded by the mean age in each 2-dimensional
bin. The bottom left panel shows the sample with PMS probability $>$90%, the
bottom right panel shows the sample wihth $>70$% probability. They are plotted
over the greyscale HR diagram of a random subset of the full evaluation
sample.

As $\rho$ Oph is very young and deeply embedded in a dusty cloud, the line of
sight extinction rises significantly behind the cloud, offering an excellent
test of the sensitivity of contamination to reddening (Figure 6). Similarly,
the cloud has a very particular shape, with several long filaments protruding
away from the center - these filaments are not actively involved in star
formation. Even without considering their distance, contamination due to
extinction can be apparent if the identified PMS candidates follow the outline
of the cloud too well.

With this in mind, we imposed several criteria to evaluate the different
trained classifier models with different hyperparameters. Ordering the model
predictions from the highest PMS probability to the lowest, for sources within
the box of $345<l<360^{\circ}$ and $15<b<25^{\circ}$, we identified the
typical probability threshold for each model where the number of sources with
$\pi<5$ begins to match the number of sources with $\pi>5$ in a given
probability bin (i.e., the point where the rate contamination/false positives
is comparable to the rate of adding bona fide PMS stars/true positives). The
best model needs to

  * •

Maximize the overall number of $\pi>5$ sources that the model identifies up to
that point

  * •

Minimize the ratio [number of sources with $\pi<5$]/[number of sources with
$\pi>5$], i.e., minimize the overall contamination fraction up to that point.

Although these criteria have been optimized for the selection of the Upper Sco
sources, they generally yield a good selection of PMS candidates in other
nearby star forming regions as well.

Selection of regions beyond 1 kpc presents a bigger challenge, as their lower
mass members tend to be too faint to be within the sensitivity limits. Thus,
members of various star forming regions located at those distances tend to
have overall lower probability than their more nearby counterparts. Because of
this, it is difficult to estimate the contamination among them.

However, almost all of the identified YSOs beyond 1 kpc should be located
close to the Galactic plane. While this is generally the case, due to the
scanning law of Gaia, there is a slight excess of (most likely false positive)
PMS candidates that are found towards the ecliptic poles in almost all models.
Such false positive are usually fainter stars that do not have 2MASS
photometry. Depending on the exact limiting threshold, this usually amounts to
a few hundred stars. Therefore, in evaluation of the best model we also
consider maximizing the number of sources with $\pi<1$ and $|b|<15$ and
minimizing the number of sources with $\pi<1$ and $|b|>20$.

With all of these considerations, we tested $>100$ different models that were
trained using different hyper-parameters and with slightly different
architectures. Most models had comparable performance in the evaluation
sample, although some had a greater difficulty in separating false positives
from true positives.

Of all of them, however, one model had almost an order of magnitude better
performance than the rest in the combined evaluation metric, although it is
unclear if it was due to the most optimal tuning of the hyperparameters, or
luck in the process of the stochastic gradient descent. Regardless, this
classifier model was chosen to be implemented into Sagitta. The HR diagram
showing the outputs of this model is shown in Figure 7.

We note that while some data driven approaches (such as, e.g., decision trees)
may be reduced to a human-readable set of conditions by which classification
takes place, this is not the case with deep learning. As such, while the model
can differentiate PMS and non-PMS stars based on their fluxes (most likely
noting in some fashion that PMS stars tend to be redder and/or over-luminous
than the main sequence stars, but not in the parameter space inhabited by red
giants), it is difficult to express precisely how these fluxes are utilized by
the model. Although beyond the scope of the current study, one could employ
model interpretability methods (e.g. Zeiler & Fergus (2014)) in an attempt to
gain some insight into model behavior.

To ensure that no residual biases in distance from massive populations
propagate to the model, we examine a synthetically generated sample of stars
based on the sample from Paper II, similar to the one described in Section
2.1. All of the synthetic stars have randomly drawn distance and extinction,
and include both the stars labeled as PMS as well as those that were more
evolved. There are some difference in the fraction of sources recovered at the
distances between 20 to 5000 pc relative to the input sample (e.g., smaller
fraction of more distant sources is recovered at the same probability
threshold, in part due to a smaller fraction of low mass stars). However, in a
uniform sample the model does not systematically favor a specific set of
distances corresponding to, e.g., the distance of the Orion Complex or Sco Cen
OB2, either in a form of better recovery of true positives than for stars at
other distances, or in form of contamination from evolved stars (Figure 8).

Figure 8: Recovery of PMS stars in a synthetically generated sample as a
function of distance. A similarly smooth distribution (with a higher recovery
fraction at a cost of more false positives) can be seen at lower thresholds as
well.

###  3.5 Age Regressor

Similarly to the classifier, the regression network to predict ages was
trained on the six photometric bands, $\pi$, and $A_{V}$. In addition to the
Gaia and 2MASS bands, we originally considered including the photometry from
the AllWISE catalog as well, but it was determined to be too noisy.

In constructing the catalog, we used 30% real data, and 70% augmented data
scattered across different distances and extinctions, for a total of a total
of $\sim 187,000$ sources. 80% of this catalog was used as a training set, 10%
was used as the development set, and 10% was withheld as a test set. The
training was done using the Adam optimizer, mean square error loss, a learning
rate of $10^{-3}$, and a batch size of 20,000 sources. Every 10 epochs we
evaluated the performance against the development set, to ensure that the
network is not overfitting on the training set (e.g., ‘memorizing’ it), but
rather that it is learning patterns that generalize to the previously unseen
data.

The training continued for $\sim$10,000 epochs. Afterwards, we continued to
train the model for $\sim$2,000 epochs on real data only, to minimize
potential artefacts that may be present in the augmented sample. However, we
note that in evaluating the ages on the test set, there were no significant
systematic differences between the models with and without the additional
2,000 epochs on real data. Similarly, in various experiments, few combinations
of hyperparameters were tested, but they tended to have comparable outputs
with few obvious differences in performance (in contrast to the experiments
with hyperparameters in the classifier). Instead, for the age regression, the
biggest gains in performance were a result of careful vetting of the labels in
the training sample.

In general, the trained model is able to qualitatively reproduce the average
ages of the populations in which stars are found (Section 4.2). It does
improve on the ability to infer the star forming histories of different
regions compared to the training sample, where usually only a single age per
population was available.

We are able to benchmark the estimates of ages for some of the stars that have
been previously observed by APOGEE spectrograph. A recent study by Olney et
al. (2020) has been able to extract calibrated $\log g$ estimates for the pre-
main sequence stars, which can be used as a proxy of age. The overall trend in
Figure 9 does show that, as expected, $\log g$ is increasing as stars evolve
and approach the main sequence.

Figure 9: Comparison of the estimated ages from Sagitta vs $\log g$ inferred
from the APOGEE spectra from Olney et al. (2020). The overdensities correspond
to specific discrete clusters targeted by APOGEE. The lines show the
theoretical PARSEC isochrones (Marigo et al., 2017) for stars with mass from
0.4 (purple) to 1 M⊙(red).

###  3.6 Uncertainties

Figure 10: Errors produced for sources with YSO probability $>$70% by
scattering model inputs by their reported uncertainties.

Convolution neural networks by default do not account for the uncertainties in
the data, nor do they output the corresponding uncertainties in the
predictions. Although some other machine learning architectures may be able to
better positioned to learn the average uncertainty in the data and provide a
resulting Bayesian posterior distribution, even they struggle to accurately
parse individual source-by-source input-by-input uncertainty, such that can be
present in e.g., photometry or in parallaxes.

As an alternative, we utilize the method used by Olney et al. (2020) and that
was used in Paper II, by generating 1000 samples per each star for the sources
in Table 4, where all the inputs are scattered by adding errors to them drawn
from a normal distribution with the width corresponding to the reported
uncertainty. Each one of these realizations of the same star was passed
through the network. Uncertainties were estimated by calculating the standard
deviation of the outputs. By using this method, our uncertainties are
indicative of both the model’s stability at any given photometric regime, and
the underlying photometric errors present in the input data. If photometric
errors were not available in any bandpass (this is occasionally common for
2MASS data, even if fluxes themselves are available), they were assigned an
uncertainty of 0.1 mag.

Despite the efficiency of neural networks, it is still time consuming to
process the entire Gaia catalog even once, let alone several times,
particularly on the machines without GPU acceleration. Thus the statistics
were generated only on the subset of the evaluation sample that has been
classified with unaltered inputs with probability of being PMS $>70$%.

When comparing the classification outputs for the unaltered sample vs the mean
classification from 1000 altered iterations from considering the
uncertainties, the mean classifications do appear to be somewhat more
accurate, better able to filter out the suspicious sources (such as those
described in Section 3.4.2 as likely false positives). The mean classification
probability is also typically lower – thus we are not likely to miss a
significant number of sources for which the mean classification returns a
higher confidence.

As expected, the scatter/errors increase for sources with fainter and more
uncertain fluxes, as well as for those that are more distant and have more
uncertain parallaxes. The scatter in the classification outputs also increases
for sources with lower certainties (as a consequence, those that are older).

The computed uncertainties in age are typically on the order of 0.1 dex
(Figure 10). They do not strongly depend on age or color of a star, although
there is a slight dependence on distance, with the nearby populations having
somewhat lower uncertainties.

##  4 Evaluation and Validation

###  4.1 Overall performance

Table 4 contains the catalog of the sources in the evaluation catalog that can
be identified as PMS sources with at least 70% confidence. This catalog
consists of 197,315 sources. Figure 11 shows the distribution of the
identified sources along the sky, according to the different cuts in
confidence levels, color-coded by their estimated ages. Figure 12 shows the
spatial distribution of stars at different age slices, to better highlight the
star forming history of the solar neighborhood.

Figure 11: The distribution stellar positions of the evaluation sample up to
PMS probabilities of 95%, 90%, and 85%. The plots are in the Galactic
coordinates and are color coded by the predicted age, and the first panel has
been annotated to indicate notable star forming regions. Note that older
sources are more apparent at lower certainty thresholds, as they are located
closer to the main sequence.  Table 4: Sample of output catalog based on DR2,
with age, classification, and predicted extinction.

Gaia DR2 | $l$ | $b$ | Predicted Age | PMS | Predicted $A_{V}$
---|---|---|---|---|---
Source ID | (deg.) | (deg.) | (dex) | Probability | (3d pos., mag.)
25220048169856 | 176.382 | -48.739 | 7.369 $\pm$ 0.177 | 0.740 $\pm$ 0.100 | 0.382 $\pm$ 0.008
2194034672517070720 | 97.309 | 9.949 | 6.817 $\pm$ 0.028 | 0.762 $\pm$ 0.131 | 0.822 $\pm$ 0.008
5224626096939825024 | 298.281 | 14.145 | 6.106 $\pm$ 0.032 | 0.955 $\pm$ 0.009 | 0.378 $\pm$ 0.001

††footnotetext: Only a portion shown here. Full table is available in an
electronic form.

Figure 12: Distribution of stars at various age ranges for sources with 85%
confidence threshold, color coded by parallax, to demonstrate the various
features that emerge for different epochs of star forming history. Note the
ring-shaped structure apparent in the last panel, discussed in section 5.2.

The evaluation catalog extends up to 5 kpc in parallax, however, at larger
distances, increasingly fewer low mass stars can be detected. Therefore, 90%
of all sources classified as PMS sources are located within 1 kpc, and 70% are
located within 500 pc. The difference becomes more extreme at particular age
ranges. Only lower mass stars can be still be identified as PMS at older ages
(e.g., $>$30 Myr), thus, the ability to identify them at larger distances is
suppressed compared to younger (e.g., $<5$ Myr) stars. Similarly, the
confidence with which PMS stars can be identified tends to be lower both for
sources that are older as well as for sources that are more distant (compared
to younger nearby sources), as both of them are located preferentially closer
to the main sequence and/or the red giant branch.

We note that only $\sim$30,000 out of $\sim$200,000 candidate PMS stars
presented here were also used in the training and testing sample from Paper
II. Furthermore, while some of them may have also been previously included in
other studies of pre-main sequence stars (most notably in Zari et al., 2018,
see Section 4.3.1), the bulk of the catalog are new identifications.

Figure 13: Left: HR diagram of the evaluation sample, color-coded by the
average probability of a source being PMS in each 2-dimensional bin. On top of
it, in blue, are the stars identified as double lined spectroscopic binaries
in the APOGEE data (M. Kounkel, et al, in prep). These sources preferentially
trace the photometric binary sequence. Right: Distribution of probabilities of
classifying these SB2s as pre-main sequence. Note, although fields associated
with notable star forming regions (e.g, Orion, Upper Sco) have been excluded
form the sample, some bona-fide PMS stars may have still persist in other
fields.

Few notes of caution should be given to the unresolved binaries. The
classifier largely avoids the binary sequence of evolved stars, particularly
at higher probabilities (Figure 13). No confusion occurs in the sources
younger than $\sim$40–60 Myr in the sources, as they are located above the
binary sequence. However some of contamination from binaries can be present at
lower probabilities (e.g., at thresholds $<$70–80%). For example, at 70%
threshold there is a hint of overdensity from the binary stars in Praesepe,
which is a 600 Myr cluster that does not contain PMS stars. At 80% threshold,
this overdensity disappears.

Table 5: Sample of output catalog based on EDR3, with age, classification, and
predicted extinction.

Gaia EDR3 | $l$ | $b$ | Predicted Age | PMS | Predicted $A_{V}$
---|---|---|---|---|---
Source ID | (deg.) | (deg.) | (dex) | Probability | (3d pos., mag.)
110951890948096 | 176.087 | -48.388 | 0.871$\pm$0.019 | 7.421$\pm$0.035 | 0.415$\pm$0.008
265532058523264 | 175.988 | -47.522 | 0.767$\pm$0.054 | 7.435$\pm$0.018 | 0.436$\pm$ 0.027
313670051618560 | 174.933 | -48.774 | 0.831$\pm$0.102 | 7.872$\pm$0.215 | 0.646$\pm$ 0.050

††footnotetext: Only a portion shown here. Full table is available in an
electronic form.

The reason for this is that the sources that are on the binary sequence have
been included in the training sample, thus the classifier has learned that the
colors that correspond to these main sequence binaries are likely false
positives. It cannot effectively separate true single PMS stars older than 60
Myr that overlap with the binary sequence. However, as such young stars are
rare in comparison to main sequence stars, the classifier downweights them
both equally to minimize the loss. Thus, minimizing contamination from main
sequence binaries results in a lack of stars $>$60 Myr included in Table 4.
Nonetheless, with independently derived membership (such as from analyzing the
distribution of stars in the phase space in older regions regions with bona
fide PMS stars, like in $\alpha$ Per), it is nonetheless possible to estimate
their ages without relying on the classifier.

Unfortunately, however, Sagitta is unable to separate pre-main sequence
unresolved binaries from single PMS stars stars. In young populations where
there is a clearly defined binary sequence for the cluster, the stars on that
binary sequence get assigned preferentially younger ages by $\sim$0.1-0.15 dex
than the age of the single stars in the same cluster. This is consistent with
the relative ages for single and binary stars that could be estimated with
traditional isochrone fitting. Thus, the ages of sources that could be
suspected to be unresolved binaries or tertiaries (such as in the cases where
binary sequence for a cluster is apparent, which can be seen in regions as
young as 8 Myr with a mono-age population e.g., Bouma et al., 2020, and
especially in the younger regions without a clearly defined binary sequence)
should be treated with care.

Recently released Gaia EDR3 (Gaia Collaboration et al., 2020) has changed the
definition of the bandpasses compared to DR2. The parallax has also been
improved by $\sim$30%. There are no strong systematics in the performance of
Sagitta when applied to EDR3, and the measurement of age and classification
compared to DR2 is generally consistent with each other within 1$\sigma$
according to the reported errors. Applying the pipeline to the full EDR3 that
meet the same quality checks as described in Section 2.3 results in a larger
catalog of stars that can be identified as likely PMS, but that is primarily
driven by these sources previously not meeting the required precision in
parallax and/or fluxes. The sources that are newly identified as PMS in EDR3
tend to be fainter and be located at preferentially larger distances,
extending the sensitivity limits of the survey. The resulting catalog is
included in the Table 5. We note that while the analysis in this paper,
including the subsequent sections, is limited to the DR2 data, overall
conclusions are consistent in the sample derived from EDR3.

###  4.2 Properties of known star forming regions

Figure 14: Distribution of stars identified by Sagitta towards various star-
forming regions, color-coded by stellar ages. Probabilities have been selected
between 85% and 95% to best represent each region, based on the age and
distance to each. To remove contamination from other nearby regions, the ages
for Taurus and Vela have been restricted to $<$ 15 Myr, and to highlight the
more distant region of Cygnus, distance was restricted to $>500$ pc.

Figure 14 shows the zoom-in view of various star forming regions that were
used as a benchmark in evaluating the measured stellar ages. Figure 15 shows
the corresponding distribution of ages for each individual SFR view.

Figure 15: Distribution of ages for star-forming regions in Figure 14. Note
that the membership of each region was not assessed beyond the position on the
sky and some distance and/or age ranges, therefore contamination from
unrelated stars is present in each field of view.

####  4.2.1 Orion Complex

The Orion Complex is the closest region of ongoing massive clustered star
formation, containing $>10,000$ stars with ages of $<1$ Myr to $>10$ Myr. The
youngest stars are found in the Orion A and B molecular clouds and older stars
are found in Orion D and $\lambda$ Ori (Kounkel et al., 2018). Sources from
Orion represent a significant fraction of the input training set, and provide
a valuable evaluation for the performance of both the classifier and the
regression model.

The model recovers $\sim$7,500 sources with $>$95% threshold, and $\sim$9,000
within $>$90% threshold. The age gradients and the typical ages that have
previously been observed by Kounkel et al. (2018) are well recovered by the
model.

We note that the evaluation catalog has a hole near Trapezium in the Orion
Nebula, as nebulosity degrades the photometric quality. Thus, few sources met
the cuts specified in Section 2.3. Applying Sagitta on a separate catalog that
has not been as constrained by the quality of inputs, it is able to recover
both the members and the appropriate ages.

####  4.2.2 $\rho$ Oph and Sco–Cen

The $\rho$ Oph star-forming region has already been discussed as a means of
evaluating the performance of the classifier, but both $\rho$ Oph and the
surrounding, slightly older Upper Scorpius region are also notable for their
peculiar star forming history.

We recover the typical age of $<1$ Myr for $\rho$ Oph, and $\sim$5 Myr for
Upper Sco (Preibisch & Mamajek, 2008; Pecaut & Mamajek, 2016). Although the
transition between the two region is rather sharp, without a significant age
spread larger than $\pm$1 Myr, there is some overlap between the two
populations, furthermore, the eastern part of Upper Sco contains somewhat
younger stars. From the Upper Sco, along the rest of Sco–Cen, we recover a
relatively smooth gradient in age 15–20 Myr towards the Lower Centaurus Crux
that has previously been observed in other works. Members of the younger Lupus
clouds (such as III and IV) get recovered with the characteristic age of
$\sim$3 Myr. The southern portion of CrA that is still associated with the
molecular gas has a typical age of $\sim$4 Myr; the northern portion that has
since cleared its gas has an age of $\sim$6–8 Myr. Nearby, there is also
$\sim$20 Myr population that also appears to be related to Sco–Cen.

While the average ages we measure for these regions are consistent to the
literature values (as the literature values were originally used for
training), we note that are able to go from discrete region-specific estimates
to a more homogeneous map of star forming history.

####  4.2.3 Vela

Towards Vela there are two unrelated populations found in a similar volume of
space. One is Vela OB2, which is associated associated with $\gamma$ Velorum,
and has a typical age of $\sim$10 Myr. The other populations has an age of
$\sim$30–35 Myr, and it contains an open cluster NGC 2547 (Jeffries et al.,
2009, 2014; Beccari et al., 2020).

Due to its youth, we recover Vela OB2 at higher confidence – the bulk of the
members can be identified at the threshold of $>$95%, containing $\sim$2,700
stars. On the other hand, NGC 2547 becomes apparent only with the threshold
lowered to $\sim$85%, containing $\sim$10,000 stars. We recover the average
ages of both these populations.

Vela OB2 in particular shows curious star forming history. The stars that are
located towards the northern group H (Cantat-Gaudin et al., 2019) are
preferentially younger than those near the central part.

####  4.2.4 Taurus and Perseus

Due their proximity and youth, the Taurus molecular clouds contain some of the
best studied young stars. This region does not contain any clusters, rather,
it is a collection of several diffuse clouds, some of which are up to 30 pc
apart. The most up-to date membership of this region is presented in Luhman
(2018). We are able to recover most of these members within our evaluation
catalog, as well as add a number of new candidates.

We are able to recover the typical age of $\sim$1–3 Myr for much of the
previously known members. There have also been a suggestion of an older nearby
$\sim$16 Myr population (Kraus et al., 2017), which we are also able to
recover. As is the case with the younger stars, this population is an assembly
of diffuse clumps of stars, resembling evolving cirrus clouds.

Along a similar sight line as Taurus (but at a somewhat larger distance) lies
Perseus. We are able to recover the age of $\sim$6–8 Myr for Per OB2 (with
some substructure), as well as ages of 1–3 Myr for younger clusters IC 348 and
NGC 1333 (Azimlu et al., 2015).

####  4.2.5 Serpens

Serpens contains several young clusters located at a distance of $\sim$450 pc.
Similarly to the work of Herczeg et al. (2019), we are able to recover stars
towards Serpens Main (with the age of $\sim 2$ Myr), as well as Serpens
Northeast, and Serpens far-South (with the age of $\sim$3–5 Myr). There also
appears to be substantial diffuse population 5–8 Myr population surrounding
them. Unfortunately, W40, likely the youngest region in this star forming
complex, is too deeply embedded to be seen in the optical regime.

We recover $\sim$500 sources towards Serpens up to the threshold of $>$95%,
and $\sim$1500 the threshold of $>$90%.

####  4.2.6 Cygnus

The star-forming regions in Cygnus, (particularly Cygnus OB2) are located at
much larger distance than other regions discussed in this work. Furthermore,
it is located behind a considerable layer of extinction (Wright et al., 2016).
Because of this, it can only be recovered in full at lower thresholds from the
classifier. Nonetheless, we are able to recover ages of $\sim$4–10 Myr for Cyg
OB2.

###  4.3 Comparison to Other Catalogs

####  4.3.1 Zari et al. (2018)

Zari et al. (2018) identified pre-main sequence stars younger than 20 Myr
within 500 pc using Gaia DR2 data. The available catalog provides three
confidence intervals for pre-main sequence sources based on their kinematic
distributions. Their total sample contains 43,719 sources, with 23,686
satisfying the strictest kinematic threshold.

Figure 16: The top panel shows a Hertzsprung-Russel Diagram constructed from
young stars in Zari et al. (2018), analyzed with Sagitta. Red points represent
sources that Sagitta assigns $<70$% likelihood of being pre-main sequence. The
bottom panel shows the spatial distribution for the aforementioned sources
with $<70$% likelihood of being young, color coded by parallax.

In examining full catalog from Zari et al. (2018), Sagitta confirms 18,488 of
the objects they identify as pre-main sequence with $>$90% probability, and
23,115 with $>$80% (we note that in the full evaluation sample there are
30,689 stars younger than 20 Myr at 90% probability within 500 pc, and 44,225
at 80% probability). There is a particularly good agreement in the identified
stars located within $\sim$200 pc within the appropriate age bounds. The
sources to which we assign lower confidence in classifications with Sagitta
($<$70%) are preferentially located at the distance of $>$400 pc (close to the
maximum distance of the Zari et al. (2018) study), and they tend to be bluer
($G_{BP}-G_{RP}<2$). Most of these sources do not appear to trace known star
forming regions, rather they trace the extinction patterns. Similarly, they do
not have much coherence in their radial velocities, as would be expected in
young populations. It is likely that these sources are contamination from the
main sequence in their sample (single or unresolved visual binary) due to
extinction (Figure 16). However, it is possible that some of them do trace
some patterns in the star forming history through massive stars (to which we
lose sensitivity) that we cannot trace with lower mass counterparts.

####  4.3.2 Marton et al. (2016, 2019)

Marton et al. (2019) used Random Forest to classify pre-main sequence sources
using Gaia DR2 and WISE photometry. Their catalog includes classifications for
101,838,724 Gaia sources, with 1,509,781 located within 5 kpc and having
greater than 90% likelihood in being classified as YSOs. The classification
was done only on the areas of the sky above a given opacity threshold using
the Planck dust opacity map. On the surface level, this catalog recovers the
underlying shape of various star forming regions (such as Orion A & B, $\rho$
Oph, Taurus-Auriga, and others). However, this is somewhat misleading, as the
opacity threshold pre-selects molecular clouds, and masks out young
populations that are no longer associated with the molecular gas. Within that
mask, however, the catalog is prone to contamination, even at a very high
level of reported certainty.

As mentioned in Section 3.4.2, $\rho$ Oph is a particularly useful region for
evaluating contamination. We compare the performance of Sagitta versus the
Marton et al. (2019) catalog in Figure 17. The left panel shows the
distribution of parallaxes of both models identified with likelihood $>95$%
towards that star forming region. Notably, while the classification from
Marton et al. (2019) does recover some sources at the appropriate parallax for
the region ($\sim$7.5 mas), the vast majority of the stars classified with
high confidence in their Gaia-ALLWISE model as YSOs have distances that are
more consistent with reddened background stars. In contrast, with the same
confidence threshold, our classification identifies a larger number of bona-
fide YSOs in the appropriate parallax range overall, with only a small degree
of contamination. Examining the Sagitta classifications using different
confidence thresholds, and using the parallax information to assess the
reliability of the classification, we see the number of likely contaminants
decrease and the number of bona fide PMS candidates increase as the threshold
increases. On the other hand, in the Gaia-ALLWISE catalog, the overall
fraction of contaminants to bona fide YSO candidates remains mostly flat
throughout the entire probability distribution, and a marginal degree of
confidence is not achieved until $>99$%.

A similar situation persists in the other star forming regions. The Gaia-
ALLWISE catalog recovers many more candidate members of these regions than
prior censuses have found, spread mostly uniformly within the outline of the
clouds, regardless of the intrinsic underlying density distribution of those
populations. As such, it is likely that their machine learning approach shows
overreliance to extinction as a proxy for youth.

In total, our evaluation catalog recovers only 11,654 sources at $>$70%
confidence threshold in common with the catalog of Marton et al. (2019) for
the sources they classify as YSOs at 90% confidence threshold.

We note that in large part, the contamination in the Gaia-ALLWISE catalog is
driven by noisy (and possibly partially mislabeled) data in the training
sample, which then propagates to the noise in the predictions (to a lesser
degree, this is also the issue in the work of Chiu et al., 2020, for distant
sources). The issue is further compounded in applying a trained model to very
uncertain data. When we analyze the full catalog of YSOs from Marton et al.
(2019) that they classify at 90% confidence, out of 1.7 million stars, Sagitta
would assign 1/4 of them $>70$% confidence, and 1/10 of them $>90$%
confidence, and the outputs from Sagitta would also be strongly contaminated
by red giants with a large parallax error. As various data-driven classifiers
are becoming more commonplace, this highlights the importance of rigorous
vetting of the data that are processed by machine learning algorithms (both in
training and in evaluating), to fully understand the limitations and ensure
that these algorithms are not applied indiscriminately.

Figure 17: Comparison of distribution of parallaxes for sources towards $\rho$
Oph and Upper Sco star-forming region as identified by Marton et al. (2019)
and Sagitta (this work) at different thresholds. Sources with small parallaxes
are likely contamination due to extinction.

Similarly, we examine the catalog from Marton et al. (2016), where they used
Support Vector Machine methods to identify YSO candidates from WISE data
alone, as that work predates the release of Gaia DR2. This catalog contains
133,980 Class I & II and 608,606 Class III YSOs, most of which are too
reddened to be detected in the optical regime, making them difficult to
evaluate. For the sources for which parallaxes are available, a significant
fraction of them do appear to be reddened supergiants. Although,
unsurprisingly, even in the subset of sources with optical emission,
contamination in Class I & II sources is somewhat less prominent than it is in
Class III, as the former tend to have peculiar colors from the protoplanetary
disks, compared to the latter which are just naked photospheres. In total we
only recover 3,722 sources from this catalog compared to our evaluation
catalog with the threshold of $>$70%, or 17,472 without any quality checks on
the data within the same threshold.

####  4.3.3 Vioque et al. (2020)

Vioque et al. (2020) used machine learning techniques to search Herbig Ae/Be
stars using Gaia, 2MASS, Sloan, IPHAS, VPHAS+, and WISE photometry. They
identify 8,470 candidate PMS stars, 693 classical main sequence Be stars, as
well as providing a list of 1,309 sources that could belong to either type
with above 50% probability.

Their selection criteria is preferentially sensitive to the stars that are
more massive than the sources we are able to identify as PMS candidates in
this work. Furthermore, based on the availability of IPHAS and VPHAS+ data,
they are restricted to only $\sim 1^{\circ}$ within the Galactic plane.
Because of this, our classifier identifies only a few sources in common with
this catalog. From their catalog we classify 3500 as PMS with Sagitta within
80% threshold. These sources tend to be very young, with an average predicted
age of approximately 5 Myr.

####  4.3.4 Kuhn et al. (2020)

Recently, Kuhn et al. (2020) have performed a data-driven selection of dusty
YSOs in the Spitzer data across the Galactic plane. As their catalog is
primarily focused on very reddened sources that do not necessarily have good
Gaia astrometry, the overlap of their selection with our evaluation catalog is
minimal, only 456 stars. Similarly, of 36,423 sources that do have optical
counterparts, regardless of data quality, our classifier would flag only 37%
of these sources as likely pre-main sequence with $>$70% confidence.
Nonetheless, the catalog does appear to be robust and complementary to our
selection, identifying preferentially younger stars and providing a more
complete selection, particularly in the more distant star forming regions,
including some that we only barely recover (such as, e.g., Sco OB1).

Nonetheless, the age estimator part of Sagitta does appear to work well on
this catalog, resulting in an average age of $\sim$4 Myr. Furthermore, it does
appear to reveal some coherent age gradients in these star forming regions.

##  5 Discussion

###  5.1 Local Bubble & Gould’s belt

Figure 18: Distribution of PMS stars (up to a confidence threshold of 80%) in
the heliocentric rectangular reference frame, color coded by age, as a
demonstration of the young stars tracing the outline of the Local Bubble. In
X-Z and Y-Z projections, the distribution of PMS stars is largely planar,
following the tilt of the Gould’s Belt (See also Figure 19).

Gould’s belt has been a long recognizable feature of the solar neighborhood,
showing the apparent tilt of star forming regions, such as Sco-Cen, Orion,
Taurus, and Perseus relative to the Galactic plane. Over the years, there have
been a number of interpretations to the causes of this tilt, whether it is
caused by a series of supernovae eruptions (Pöppel & Marronetti, 2000), or a
collision of some sort with the disk (Comeron & Torra, 1994; Bekki, 2009).

Recently, some of the populations on one side of the Gould’s belt (such as
Orion) have been associated with the Radcliffe Wave (Alves et al., 2020) \- a
2.7 kpc long structure that has a number of ripples protruding from the
Galactic plane, up to 160 pc in amplitude. However, no interpretation has been
given as a cause for these ripples.

Zari et al. (2018) have analyzed the distribution of stars younger than 20 Myr
and found no evidence for a fully connected Gould’s belt, rather that all of
the individual populations appeared to be unrelated. Similarly Bouy & Alves
(2015) have used a census of nearby OB stars and arrived to similar
conclusions. With an improved census of PMS stars that extend towards the
older ages we seek reexamine this.

Figure 18 shows the top down map of PMS stars. There is a clear ring-like
structure with the radius of $\sim$100–150 pc that connects Sco–Cen and Taurus
as well as some of the older populations, such as $\alpha$ Per/Cas Tau OB
association. Although hints of a complete ring are seen in stars $<20$ Myr old
(particularly tracing the inner rim), it is the older stars that define this
ring most clearly, which may be part of the reason why Zari et al. (2018) have
not identified it in their data.

Furthermore, beyond this ring, there appears to be a gap in the 3d
distribution of PMS stars at the distances of $\sim$200–250 pc, with other
populations becoming more prominent at distances of $>$300 pc.

This structure does not appear to be artificial, persisting even at higher
confidence levels. Furthermore, as discussed in Section 3.4.2, it cannot be
attributed to the classifier favoring a particular set of distances in a truly
uniform distribution of stars as the classifier does not intrinsically favor
any specific distance in either the recovery fraction or in contamination
(Figure 8.)

Outside of very low mass populations such as TW Hya (which we do recover), a
lack of significant young star forming regions within 100 pc has long since
been known. Indeed, the Sun is located near the center of the Local Bubble, a
cavity of significantly lower density of neutral hydrogen compared to what is
typically found in the interstellar medium. Using various tracers, such as
diffuse interstellar bands (that can trace the X-ray dissociation within the
Bubble Farhang et al., 2019) or extinction (Lallement et al., 2018; Leike &
Enßlin, 2019) to trace the morphology of the local bubble, it fits well within
the identified ring of stars.

It has been noted Ophiuchus and Taurus have velocities that are comparable in
magnitude, but opposite in direction, and they would have been in proximity of
one another $\sim$20-25 Myr ago (Rivera et al., 2015). This trace back age is
comparable to the average age of the stars in the ring, although there are
also a number of $\sim$40 Myr old stars that compose it. Unfortunately, so far
only a few sources have radial velocities available to confirm a global
expansion along other sight lines, or estimate the expansion age. At a glance
this appears to be an effect of a supernova explosion. It should be noted that
different populations do appear to have a somewhat different peculiar velocity
relative to one another, at least in the proper motion space. Thus most
likely, instead of shockwave clearing a gas of a particular population (as has
been the case of supernovae in young star forming regions such as Orion, and
potentially Vela; Kounkel, 2020; Großschedl et al., 2020; Cantat-Gaudin et
al., 2019) the shock front associated with the expanding bubble may have
rammed into the neighboring clouds, not dissimilar to a scenario described by
Inutsuka et al. (2015) in which molecular clouds trace interaction regions
between even shorter lived bubbles.

The Local Bubble does not immediately explain the lack of star forming regions
at the distances of $\sim$200–250 pc. Such a gap can be seen in the
distribution of present day molecular clouds (Zucker et al., 2020), and it is
also present in the catalog from Zari et al. (2018). Although it is not
impossible, it would be surprising for two unrelated events to occur in the
vicinity of the Sun to form two separate rings at $\sim 100$ pc and at
$\sim$300 pc. Instead, it may be possible that they are produced by a related
event, possibly caused by multiple shock fronts. This may suggest a common
origin for the populations that are a part of the Gould’s belt.

###  5.2 30 Myr Bubble

Figure 19: Top: 3d distribution of stars of the 30 Myr bubble relative to the
other stars in the catalog. Middle: distribution of the sources towards the
bubble in the plane of the sky, with proper motions, in the local standard of
rest, corrected for the average velocity of the stars, with the thick part of
the arrow showing the current position of the sources. Arrows are color coded
by the cosine of angle proper motions have relative to the center of the
bubble, with red (-1) showing sources moving radially inward and blue (+1)
showing those moving radially outward. Bottom: a histogram showing the
distribution of the aforementioned radial component of motion.

There is another peculiar bubble-like structure that can be identified in the
data. It can be seen in the last panel of Figure 12, primarily traced by the
stars older than 25 Myr, towards the direction of the Galactic center, and it
is $\sim 40^{\circ}$ in radius. This bubble appears to begin at the distance
of $\sim$200 pc and forming a hemisphere $\sim$400 pc in diameter (Figure 19,
top). Unfortunately the distant edge of this bubble is difficult to detect as
at those distances we lose sensitivity to low mass stars that can be used as
tracers for this age range.

This bubble is not caused by extinction towards the Galactic center. Although
there are a number of optically thick clouds in the volume of space associated
with it (e.g., the Aquilla rift, for which we do recover a sizable population
of stars in the 2–10 Myr range), the edge of the bubble is located far outside
of those clouds.

In analyzing the proper motions of the stars located on the other edge of the
bubble we identify a somewhat peculiar pattern: corrected for the average
velocity of all the stars in the sample, there is a strong preference for them
to move either directly radially inwards towards the center (at $l\sim
6^{\circ}$, $b\sim 6^{\circ}$) or radially outward away from it, with next to
no tangential component in the proper motions (Figure 19). It is unclear what
could be the cause of such a signature.

Based on the velocities of stars that are just expanding, we would estimate
the expansion age of $\sim$13 Myr, or approximately half of their age.

##  6 Conclusions

One of the outstanding questions in the star formation community is how can
post T Tauri stars be identified. We present an automated method of
identifying PMS stars and estimating their ages (up to $\sim$70 Myr) through
Gaia DR2 and 2MASS photometry using a neural network framework. This allows
for a homogeneous analysis of large volumes of data characterizing star
forming history of the solar neighborhood. Furthermore, this approach is not
reliant on a kinematic selection, making it possible to search for
kinematically peculiar young stars, such as runaways (e.g., McBride & Kounkel,
2019).

Applying a classifier to a curated subset of Gaia DR2 data with most reliable
astrometry and photometry, we identify 197,315 stars as likely PMS sources
with confidence of $>70$%, and 496,488 stars in Gaia EDR3 data. The code is
made available on GitHub222https://github.com/hutchresearch/Sagitta, to enable
the usage of Sagitta outside of this curated subset.

Sagitta is robust against contamination, especially when compared a number of
previous studies that also aimed to identify young stars using optical and
near infrared data. The precise confidence threshold that should be used to
select candidate PMS stars in a particular region depends on the distance to
and the age of the population that one seeks to characterize, however.

The estimated ages that we measure are consistent with what has been
previously measured in some of the better studied star forming regions.
Furthermore, in many cases, they allow for a more granular look at the
evolution of various populations than what was previously available. It should
be noted, however, that caution should be expressed regarding the pre-main
sequence binaries, as they may appear to be systematically younger than they
are.

In examining the distribution of stars in the solar neighborhood, we identify
various features. Most notably we identify a ring of stars at $\sim$100 pc
with ages of up to $\sim$40 Myr, tracing the outer edges of the Local Bubble.
It is likely that the formation of this bubble have lead to the formation of
the Gould’s belt. We also find a second bubble consisting of $\sim$30 Myr old
stars in the direction towards the Galactic center.

In future, a follow up of the sample presented in this work by large
spectroscopic surveys (such as SDSS-V APOGEE) would be of great benefit to
confirming the candidates, as well as allowing for a better understanding of
the dynamical and chemical evolution of PMS stars.

We thank Eleonora Zari and Keivan Stassun for giving feedback on the
manuscript. A.M, M.K., and K.C. acknowledge support provided by the NSF
through grant AST-1449476. This work has made use of data from the European
Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed
by the Gaia Data Processing and Analysis Consortium (DPAC,
https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has
been provided by national institutions, in particular the institutions
participating in the Gaia Multilateral Agreement. The authors thank the Nvidia
Corporation for their donation of GPUs used in this work.

## References

  * Alves et al. (2020) Alves, J., Zucker, C., Goodman, A. A., et al. 2020, Nature, doi: 10.1038/s41586-019-1874-z
  * Andrews & Wolk (2008) Andrews, S. M., & Wolk, S. J. 2008, Bo Reipurth, I, 1.  https://arxiv.org/abs/arXiv:0809.2276v1
  * Azimlu et al. (2015) Azimlu, M., Martínez-Galarza, J. R., & Muench, A. A. 2015, AJ, 150, 95, doi: 10.1088/0004-6256/150/3/95
  * Bally (2008) Bally, J. 2008, Overview of the Orion Complex, ed. B. Reipurth, 459
  * Baraffe et al. (1998) Baraffe, I., Chabrier, G., Allard, F., & Hauschildt, P. H. 1998, A&A, 337, 403
  * Baraffe et al. (2015) Baraffe, I., Homeier, D., Allard, F., & Chabrier, G. 2015, A&A, 577, A42, doi: 10.1051/0004-6361/201425481
  * Beccari et al. (2020) Beccari, G., Boffin, H. M. J., & Jerabkova, T. 2020, MNRAS, 491, 2205, doi: 10.1093/mnras/stz3195
  * Bekki (2009) Bekki, K. 2009, MNRAS, 398, L36, doi: 10.1111/j.1745-3933.2009.00702.x
  * Bouma et al. (2020) Bouma, L. G., Winn, J. N., Ricker, G. R., et al. 2020, arXiv e-prints, arXiv:2005.10253.  https://arxiv.org/abs/2005.10253
  * Bouy & Alves (2015) Bouy, H., & Alves, J. 2015, A&A, 584, A26, doi: 10.1051/0004-6361/201527058
  * Briceño et al. (2019) Briceño, C., Calvet, N., Hernández, J., et al. 2019, AJ, 157, 85, doi: 10.3847/1538-3881/aaf79b
  * Cantat-Gaudin et al. (2019) Cantat-Gaudin, T., Mapelli, M., Balaguer-Núñez, L., et al. 2019, A&A, 621, A115, doi: 10.1051/0004-6361/201834003
  * Carpenter & Hodapp (2008) Carpenter, J., & Hodapp, K. 2008, I, 1.  https://arxiv.org/abs/0809.1396
  * Chiu et al. (2020) Chiu, Y.-L., Ho, C.-T., Wang, D.-W., & Lai, S.-P. 2020, arXiv e-prints, arXiv:2007.06235.  https://arxiv.org/abs/2007.06235
  * Choi et al. (2016) Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102, doi: 10.3847/0004-637X/823/2/102
  * Cohen & Kuhi (1979) Cohen, M., & Kuhi, L. V. 1979, ApJS, 41, 743, doi: 10.1086/190641
  * Comerón (2008) Comerón, F. 2008, Handbook of Star Forming Regions, I, 295.  http://www.aspmonographs.org/a/volumes/article_details/?paper_id=43
  * Comeron & Torra (1994) Comeron, F., & Torra, J. 1994, A&A, 281, 35
  * Covey et al. (2010) Covey, K. R., Lada, C. J., Román-Zúñiga, C., et al. 2010, ApJ, 722, 971, doi: 10.1088/0004-637X/722/2/971
  * Da Rio et al. (2012) Da Rio, N., Robberto, M., Hillenbrand, L. A., Henning, T., & Stassun, K. G. 2012, ApJ, 748, 14, doi: 10.1088/0004-637X/748/1/14
  * Dahm (2008) Dahm, S. E. 2008, The Young Cluster and Star Forming Region NGC 2264, ed. B. Reipurth, 966
  * Damiani et al. (2019) Damiani, F., Prisinzano, L., Pillitteri, I., Micela, G., & Sciortino, S. 2019, A&A, 623, A112, doi: 10.1051/0004-6361/201833994
  * D’Antona & Mazzitelli (1994) D’Antona, F., & Mazzitelli, I. 1994, ApJS, 90, 467, doi: 10.1086/191867
  * Delgado et al. (2011) Delgado, A. J., Alfaro, E. J., & Yun, J. L. 2011, A&A, 531, A141, doi: 10.1051/0004-6361/201116491
  * Erickson et al. (2015) Erickson, K. L., Wilking, B. A., Meyer, M. R., et al. 2015, AJ, 149, 103, doi: 10.1088/0004-6256/149/3/103
  * Evans et al. (2009) Evans, II, N. J., Dunham, M. M., Jørgensen, J. K., et al. 2009, ApJS, 181, 321, doi: 10.1088/0067-0049/181/2/321
  * Fang et al. (2013) Fang, M., Kim, J. S., van Boekel, R., et al. 2013, ApJS, 207, 5, doi: 10.1088/0067-0049/207/1/5
  * Fang et al. (2017) Fang, M., Kim, J. S., Pascucci, I., et al. 2017, AJ, 153, 188, doi: 10.3847/1538-3881/aa647b
  * Farhang et al. (2019) Farhang, A., van Loon, J. T., Khosroshahi, H. G., Javadi, A., & Bailey, M. 2019, Nature Astronomy, 3, 922, doi: 10.1038/s41550-019-0814-z
  * Farias et al. (2020) Farias, J. P., Tan, J. C., & Eyer, L. 2020, ApJ, 900, 14, doi: 10.3847/1538-4357/aba699
  * Fischer et al. (2017) Fischer, W. J., Megeath, S. T., Furlan, E., et al. 2017, ApJ, 840, 69, doi: 10.3847/1538-4357/aa6d69
  * Gaia Collaboration et al. (2020) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2020, arXiv e-prints, arXiv:2012.01533
  * Gaia Collaboration et al. (2018) —. 2018, A&A, 616, A1, doi: 10.1051/0004-6361/201833051
  * Galli et al. (2019) Galli, P. A. B., Loinard, L., Bouy, H., et al. 2019, A&A, 630, A137, doi: 10.1051/0004-6361/201935928
  * Getman et al. (2014) Getman, K. V., Feigelson, E. D., Kuhn, M. A., et al. 2014, ApJ, 787, 108, doi: 10.1088/0004-637X/787/2/108
  * Green et al. (2019) Green, G. M., Schlafly, E., Zucker, C., Speagle, J. S., & Finkbeiner, D. 2019, ApJ, 887, 93, doi: 10.3847/1538-4357/ab5362
  * Greene & Meyer (1995) Greene, T. P., & Meyer, M. R. 1995, ApJ, 450, 233, doi: 10.1086/176134
  * Großschedl et al. (2020) Großschedl, J. E., Alves, J., & Meingast, S. 2020, arXiv e-prints, arXiv:2007.07254.  https://arxiv.org/abs/2007.07254
  * Hartmann (2009) Hartmann, L. 2009, Accretion Processes in Star Formation: Second Edition (Cambridge University Press)
  * Herbig & Reipurth (2008) Herbig, G., & Reipurth, B. 2008, Handbook of Star Forming Regions, Volume I, II, 108
  * Herczeg & Hillenbrand (2014) Herczeg, G. J., & Hillenbrand, L. A. 2014, ApJ, 786, 97, doi: 10.1088/0004-637X/786/2/97
  * Herczeg & Hillenbrand (2015) —. 2015, ApJ, 808, 23, doi: 10.1088/0004-637X/808/1/23
  * Herczeg et al. (2019) Herczeg, G. J., Kuhn, M. A., Zhou, X., et al. 2019, ApJ, 878, 111, doi: 10.3847/1538-4357/ab1d67
  * Hillenbrand et al. (2008) Hillenbrand, L. A., Bauermeister, A., & White, R. J. 2008, in Astronomical Society of the Pacific Conference Series, Vol. 384, 14th Cambridge Workshop on Cool Stars, Stellar Systems, and the Sun, ed. G. van Belle, 200.  https://arxiv.org/abs/astro-ph/0703642
  * Inutsuka et al. (2015) Inutsuka, S.-i., Inoue, T., Iwasaki, K., & Hosokawa, T. 2015, A&A, 580, A49, doi: 10.1051/0004-6361/201425584
  * Jackson et al. (2018) Jackson, R. J., Deliyannis, C. P., & Jeffries, R. D. 2018, MNRAS, 476, 3245, doi: 10.1093/mnras/sty374
  * Jeffries et al. (2009) Jeffries, R. D., Naylor, T., Walter, F. M., Pozzo, M. P., & Devey, C. R. 2009, MNRAS, 393, 538, doi: 10.1111/j.1365-2966.2008.14162.x
  * Jeffries et al. (2014) Jeffries, R. D., Jackson, R. J., Cottaar, M., et al. 2014, A&A, 563, A94, doi: 10.1051/0004-6361/201323288
  * Kenyon et al. (2008) Kenyon, S. J., Gómez, M., & Whitney, B. A. 2008, Low Mass Star Formation in the Taurus-Auriga Clouds, ed. B. Reipurth, 405
  * Koenig et al. (2012) Koenig, X. P., Leisawitz, D. T., Benford, D. J., et al. 2012, ApJ, 744, 130, doi: 10.1088/0004-637X/744/2/130
  * Kos et al. (2019) Kos, J., Bland-Hawthorn, J., Asplund, M., et al. 2019, A&A, 631, A166, doi: 10.1051/0004-6361/201834710
  * Kounkel (2020) Kounkel, M. 2020, arXiv e-prints, arXiv:2007.09160.  https://arxiv.org/abs/2007.09160
  * Kounkel & Covey (2019) Kounkel, M., & Covey, K. 2019, AJ, 158, 122, doi: 10.3847/1538-3881/ab339a
  * Kounkel et al. (2020) Kounkel, M., Covey, K., & Stassun, K. G. 2020, arXiv e-prints, arXiv:2004.07261.  https://arxiv.org/abs/2004.07261
  * Kounkel et al. (2018) Kounkel, M., Covey, K., Suárez, G., et al. 2018, AJ, 156, 84, doi: 10.3847/1538-3881/aad1f1
  * Kraus et al. (2017) Kraus, A. L., Herczeg, G. J., Rizzuto, A. C., et al. 2017, ApJ, 838, 150, doi: 10.3847/1538-4357/aa62a0
  * Ksoll et al. (2018) Ksoll, V. F., Gouliermis, D. A., Klessen, R. S., et al. 2018, MNRAS, 479, 2389, doi: 10.1093/mnras/sty1317
  * Kuhn et al. (2020) Kuhn, M. A., de Souza, R. S., Krone-Martins, A., et al. 2020, arXiv e-prints, arXiv:2011.12961.  https://arxiv.org/abs/2011.12961
  * Kumar et al. (2014) Kumar, B., Sharma, S., Manfroid, J., et al. 2014, A&A, 567, A109, doi: 10.1051/0004-6361/201323027
  * Kun et al. (2009) Kun, M., Balog, Z., Kenyon, S. J., Mamajek, E. E., & Gutermuth, R. A. 2009, ApJS, 185, 451, doi: 10.1088/0067-0049/185/2/451
  * Kun et al. (2008) Kun, M., Balog, Z., Mizuno, N., et al. 2008, MNRAS, 391, 84, doi: 10.1111/j.1365-2966.2008.13898.x
  * Lallement et al. (2018) Lallement, R., Capitanio, L., Ruiz-Dern, L., et al. 2018, A&A, 616, A132, doi: 10.1051/0004-6361/201832832
  * Leike & Enßlin (2019) Leike, R. H., & Enßlin, T. A. 2019, A&A, 631, A32, doi: 10.1051/0004-6361/201935093
  * Lindegren et al. (2018) Lindegren, L., Hernández, J., Bombrun, A., et al. 2018, A&A, 616, A2, doi: 10.1051/0004-6361/201832727
  * López Martí et al. (2013) López Martí, B., Jiménez-Esteban, F., Bayo, A., et al. 2013, A&A, 556, A144, doi: 10.1051/0004-6361/201321217
  * Luhman (2008) Luhman, K. L. 2008, Chamaeleon, ed. B. Reipurth, Vol. 5, 169
  * Luhman (2018) —. 2018, AJ, 156, 271, doi: 10.3847/1538-3881/aae831
  * Marigo et al. (2017) Marigo, P., Girardi, L., Bressan, A., et al. 2017, ApJ, 835, 77, doi: 10.3847/1538-4357/835/1/77
  * Marton et al. (2016) Marton, G., Tóth, L. V., Paladini, R., et al. 2016, MNRAS, 458, 3479, doi: 10.1093/mnras/stw398
  * Marton et al. (2019) Marton, G., Ábrahám, P., Szegedi-Elek, E., et al. 2019, MNRAS, doi: 10.1093/mnras/stz1301
  * McBride & Kounkel (2019) McBride, A., & Kounkel, M. 2019, ApJ, 884, 6, doi: 10.3847/1538-4357/ab3df9
  * Megeath et al. (2012) Megeath, S. T., Gutermuth, R., Muzerolle, J., et al. 2012, AJ, 144, 192, doi: 10.1088/0004-6256/144/6/192
  * Neuhäuser & Forbrich (2008) Neuhäuser, R., & Forbrich, J. 2008, II, 1.  https://arxiv.org/abs/0808.3374
  * Olney et al. (2020) Olney, R., Kounkel, M., Schillinger, C., et al. 2020, AJ, 159, 182, doi: 10.3847/1538-3881/ab7a97
  * Palla & Stahler (2002) Palla, F., & Stahler, S. W. 2002, ApJ, 581, 1194, doi: 10.1086/344293
  * Panwar et al. (2018) Panwar, N., Pandey, A. K., Samal, M. R., et al. 2018, AJ, 155, 44, doi: 10.3847/1538-3881/aa9f1b
  * Paszke et al. (2017) Paszke, A., Gross, S., Chintala, S., et al. 2017, in NIPS-W
  * Pecaut & Mamajek (2016) Pecaut, M. J., & Mamajek, E. E. 2016, MNRAS, 461, 794, doi: 10.1093/mnras/stw1300
  * Pöppel & Marronetti (2000) Pöppel, W. G. L., & Marronetti, P. 2000, A&A, 358, 299
  * Prato et al. (2008) Prato, L., Rice, E. L., & Dame, T. M. 2008, Where are all the Young Stars in Aquila?, ed. B. Reipurth, Vol. 4, 18
  * Preibisch & Mamajek (2008) Preibisch, T., & Mamajek, E. 2008, The Nearest OB Association: Scorpius-Centaurus (Sco OB2), ed. B. Reipurth, Vol. 5, 235
  * Prisinzano et al. (2018) Prisinzano, L., Damiani, F., Guarcello, M. G., et al. 2018, A&A, 617, A63, doi: 10.1051/0004-6361/201833172
  * Prusti et al. (1992) Prusti, T., Adorf, H. M., & Meurs, E. J. A. 1992, A&A, 261, 685
  * Rauw & De Becker (2008) Rauw, G., & De Becker, M. 2008, II, 1.  https://arxiv.org/abs/0808.3887
  * Reipurth (2008a) Reipurth, B. 2008a, The Young Cluster NGC 6604 and the Serpens OB2 Association, ed. B. Reipurth, Vol. 5, 590
  * Reipurth (2008b) —. 2008b, Handbook of Star Forming Regions, Volume I: The Northern Sky, ASP Monographs No. v. 1 (Astronomical Society of the Pacific).  http://aspmonographs.org/a/volumes/table_of_contents?book_id=1
  * Reipurth (2008c) —. 2008c, Handbook of Star Forming Regions, Volume II: The Southern Sky, ASP Monographs No. v. 2 (Astronomical Society of the Pacific).  http://aspmonographs.org/a/volumes/table_of_contents?book_id=2
  * Reipurth & Schneider (2008) Reipurth, B., & Schneider, N. 2008, Handbook of Star Forming Regions, Volume I, I, 36
  * Reipurth & Yan (2008) Reipurth, B., & Yan, C. H. 2008, Star Formation and Molecular Clouds towards the Galactic Anti-Center, ed. B. Reipurth, Vol. 4, 869
  * Rivera et al. (2015) Rivera, J. L., Loinard, L., Dzib, S. A., et al. 2015, ApJ, 807, 119, doi: 10.1088/0004-637X/807/2/119
  * Román-Zúñiga & Lada (2008) Román-Zúñiga, C. G., & Lada, E. A. 2008, I, 1.  https://arxiv.org/abs/0810.0931
  * Schoettler et al. (2020) Schoettler, C., de Bruijne, J., Vaher, E., & Parker, R. J. 2020, MNRAS, 495, 3104, doi: 10.1093/mnras/staa1228
  * Siess et al. (2000) Siess, L., Dufour, E., & Forestini, M. 2000, A&A, 358, 593
  * Smith & Brooks (2008) Smith, N., & Brooks, K. J. 2008, II, 1.  https://arxiv.org/abs/0809.5081
  * Suárez et al. (2017) Suárez, G., Downes, J. J., Román-Zúñiga, C., et al. 2017, AJ, 154, 14, doi: 10.3847/1538-3881/aa733a
  * Taylor (2005) Taylor, M. B. 2005, in Astronomical Society of the Pacific Conference Series, Vol. 347, Astronomical Data Analysis Software and Systems XIV, ed. P. Shopbell, M. Britton, & R. Ebert, 29
  * Tóth et al. (2014) Tóth, L. V., Marton, G., Zahorecz, S., et al. 2014, PASJ, 66, 17, doi: 10.1093/pasj/pst017
  * Tothill et al. (2008) Tothill, N. F. H., Gagné, M., Stecklum, B., & Kenworthy, M. A. 2008, II, 1\.  https://arxiv.org/abs/0809.3380
  * Vioque et al. (2020) Vioque, M., Oudmaijer, R. D., Schreiner, M., et al. 2020, arXiv e-prints, arXiv:2005.01727.  https://arxiv.org/abs/2005.01727
  * Walawender et al. (2008) Walawender, J., Bally, J., Di Francesco, J., Jørgensen, J., & Getman, K. 2008, Bo Reipurth, I, 1.  http://www.ifa.hawaii.edu/publications/preprints/08preprints/Walawender{_}08-206.pdf
  * Wilking et al. (2008) Wilking, B., Gagne, M., & Allen, L. 2008, II, 1.  https://arxiv.org/abs/0811.0005
  * Wolk et al. (2008) Wolk, S. J., Bourke, T. L., & Vigil, M. 2008, Handbook of Star Forming Regions, II, 124.  https://arxiv.org/abs/arXiv:0808.3385v1
  * Wright et al. (2016) Wright, N. J., Bouy, H., Drew, J. E., et al. 2016, MNRAS, 460, 2593, doi: 10.1093/mnras/stw1148
  * Zari et al. (2019) Zari, E., Brown, A. G. A., & de Zeeuw, P. T. 2019, A&A, 628, A123, doi: 10.1051/0004-6361/201935781
  * Zari et al. (2018) Zari, E., Hashemi, H., Brown, A. G. A., Jardine, K., & de Zeeuw, P. T. 2018, A&A, 620, A172, doi: 10.1051/0004-6361/201834150
  * Zeiler & Fergus (2014) Zeiler, M. D., & Fergus, R. 2014, in Proc. ECCV, 818–833, doi: 10.1007/978-3-319-10590-1_53
  * Zucker et al. (2020) Zucker, C., Speagle, J. S., Schlafly, E. F., et al. 2020, A&A, 633, A51, doi: 10.1051/0004-6361/201936145