# Stars with Photometrically Young Gaia Luminosities Around the Solar System
(SPYGLASS) I: Mapping Young Stellar Structures and their Star Formation
Histories

Ronan Kerr Department of Astronomy, University of Texas at Austin
2515 Speedway, Stop C1400
Austin, Texas, USA 78712-1205
Aaron C. Rizzuto  51 Pegasi b Fellow  Department of Astronomy, University of
Texas at Austin
2515 Speedway, Stop C1400
Austin, Texas, USA 78712-1205
Adam L. Kraus Department of Astronomy, University of Texas at Austin
2515 Speedway, Stop C1400
Austin, Texas, USA 78712-1205
Stella S. R. Offner  Department of Astronomy, University of Texas at Austin
2515 Speedway, Stop C1400
Austin, Texas, USA 78712-1205

(Accepted May 16, 2021)

###### Abstract

Young stellar associations hold a star formation record that can persist for
millions of years, revealing the progression of star formation long after the
dispersal of the natal cloud. To identify nearby young stellar populations
that trace this progression, we have designed a comprehensive framework for
the identification of young stars, and use it to identify $\sim$3$\times
10^{4}$ candidate young stars within a distance of 333 pc using Gaia DR2.
Applying the HDBSCAN clustering algorithm to this sample, we identify 27 top-
level groups, nearly half of which have little to no presence in previous
literature. Ten of these groups have visible substructure, including notable
young associations such as Orion, Perseus, Taurus, and Sco-Cen. We provide a
complete subclustering analysis on all groups with substructure, using age
estimates to reveal each region’s star formation history. The patterns we
reveal include an apparent star formation origin for Sco-Cen along a
semicircular arc, as well as clear evidence for sequential star formation
moving away from that arc with a propagation speed of $\sim$4 km s-1 ($\sim$4
pc Myr-1). We also identify earlier bursts of star formation in Perseus and
Taurus that predate current, kinematically identical active star-forming
events, suggesting that the mechanisms that collect gas can spark multiple
generations of star formation, punctuated by gas dispersal and cloud regrowth.
The large spatial scales and long temporal scales on which we observe star
formation offer a bridge between the processes within individual molecular
clouds and the broad forces guiding star formation at galactic scales.

young stars — open clusters and associations — star formation — T Tauri stars
— Sco-Cen — Taurus — Perseus — Orion

††facilities: Gaia††software: Astropy (Astropy Collaboration et al., 2013),
HDBSCAN (McInnes et al., 2017), Matplotlib (Hunter, 2007), SciPy (Virtanen et
al., 2020)

##  1 Introduction

Most young stars are not found in isolation, instead residing in co-moving
star clusters or associations (e.g., Lada & Lada 2003; Krumholz et al. 2019).
These stellar overdensities are direct remnants of the molecular cloud that
collapsed to create them, and as such they preserve significant information on
the structure and dynamics of those parent clouds (Elmegreen & Lada, 1977;
Krause et al., 2020). While studies of sites of active star formation are
popular due to the presence of both young stars and dense gas for which
exquisite dynamical studies can be performed (e.g., Palla & Stahler, 1999;
Hatchell et al., 2005; Tobin et al., 2009; Kirk et al., 2013, 2017; Kerr et
al., 2019), studies of clusters and associations provide a much longer-term
view of star formation that has unique value. Rather than providing a snapshot
in time, a detailed study of a young association can trace back tens of
millions of years of star formation, enabling a complete study of a star-
forming event from its onset to present day, and by extension the processes
that drove the evolution of the population (e.g.,de Zeeuw et al. 1999; Pecaut
& Mamajek 2016; Wright & Mamajek 2018). Furthermore, as these associations
persist long after gas dispersal, the study of stellar populations through
associations may also enable the identification of more rapid forms of star
formation that do not have known equivalents in active star-forming sites.

Some processes in star formation require a complete record of the star-forming
event, not just a snapshot, to be properly studied in nature. One example is
sequential star formation, a process in which previous generations of stars
compress the cloud beside them, which can then collapse to form stars,
producing a self-sustaining cycle of star formation that can slowly propagate
across an entire molecular cloud (Elmegreen & Lada, 1977). Most cases where a
sequential process has been suggested include just two generations of star
formation: one recently-formed generation powering an H II region, and one
site of active star formation triggered in a shell that the previous
generation compressed (e.g., Lee et al., 2005; Maaskant et al., 2011; Nony et
al., 2021). Given that these processes are capable of continuing without limit
as long as unused gas remains, the current view of sequential star formation
has yet to explore large scales in both time and space.

One of the greatest strengths in using clusters and associations as a record
of past star formation lies in the time over which detailed information can be
extracted. Having such an extensive star formation record allows for the study
of long-lived star-forming processes, while also revealing unexpected
anomalies in time, such as periods of dormancy in the star formation record.
As simulations become increasingly sophisticated and begin to include more
physical processes (e.g., Grudić et al., 2020a), having a robust record of
star formation, complete with currently unexplained features can provide
critical comparisons capable of testing new theoretical frameworks.

While more spatially compact open clusters have long ago been discovered and
catalogued (e.g., Trumpler, 1921; Klein Wassink, 1927; Mermilliod, 1995), much
sparser stellar associations are considerably more challenging to identify due
to relatively sparse on-sky densities and large spatial extents, often barring
identification from on-sky density. For the nearest populations, stars with
strong indicators of youth such as protoplanetary disks have typically been
used as signposts for the identification of associated populations nearby,
such as TW Hydrae and $\beta$ Pictoris (e.g., Kastner et al., 1997; Zuckerman
et al., 2001). Beyond about 50-100 pc, however, recognizing stars with youth
indicators and finding potential companions becomes increasingly difficult,
leaving the population of low-mass, more distant associations largely
unexplored. Due to their unbound nature, small velocities relative to a
typical field star, and geometric effects from their wide spatial
distributions, the kinematics of associations are often difficult to
disentangle from the field. The effective identification of young stellar
populations therefore requires both the suppression of older background
populations and the use of accurate kinematics and 3-d spatial coordinates to
properly group the young stars together and distinguish a young association
from the field. Improvements in measurements of parallax and proper motions
can therefore significantly expand our knowledge of these associations, and
such developments frequently result in the addition of new members and the
revision of associations’ known spatial extents (e.g., de Zeeuw et al., 1999;
Preibisch & Mamajek, 2008; Rizzuto et al., 2015; Zari et al., 2018; Kraus et
al., 2017; Luhman, 2018). With its unprecedented breadth of spatial,
kinematic, and photometric data, Gaia Data Release 2 provides a data set
capable of dramatically expanding our knowledge of young stellar populations
(Gaia Collaboration et al., 2018a). Gaia photometry can be combined with
parallaxes to identify stars based on their high locations on the HR diagram,
and those stars can then be clustered according to their spatial coordinates
and transverse motions.

There have been multiple recent searches for stellar populations in the solar
neighborhood enabled by Gaia DR2, including Sim et al. (2019), which reports
the discovery of 207 new open clusters, and Kounkel & Covey (2019), which
notes the identification of 1901 stellar overdensities in Gaia DR2. However,
these surveys included stars of all ages in their clustering, so young
populations may not always stand out above the field density of the older
stars. The work of Zari et al. (2018) does focus on these younger populations,
separating out pre-main sequence stars quite effectively and revealing
multiple associations in the form of stellar overdensities. While Zari et al.
(2018) revealed substantial young stellar populations, it did not include a
clustering analysis, and as such the structures present within those
populations have yet to be rigorously defined. Some recent investigations have
focused on young stellar populations while also performing a complete
clustering analysis, however these have all been restricted to individual
associations, such as Zari et al. (2019), which studies Orion, and Cantat-
Gaudin et al. (2019b), which focuses on Vela. No paper to date has performed a
spatially unbiased, all-sky search for young stellar populations with both
robust youth diagnostics and a broad clustering analysis.

In this paper, we present the deepest comprehensive study of young stellar
populations in the solar neighborhood to date. In Section 2, we present the
Gaia DR2 data set used for this work. We outline our methods for the
identification of young stars, age estimation, and clustering in Section 3,
while a detailed analysis into the populations we identify is performed in
Section 4. Possible implications of this work with respect to the study of
star formation progression in the solar neighborhood are discussed in Section
5, while we conclude in Section 6.

##  2 Data

###  2.1 Gaia Astrometry and Photometry

In this paper, we search for nearby young stars among a large sample drawn
from Gaia Data Release 2 (Gaia Collaboration et al., 2018a; Lindegren et al.,
2018). When querying sources from the Gaia Archive, several search
restrictions are included to ensure a manageable sample size, and to guarantee
that all stars have high-quality photometric and astrometric measurements.
First, we impose a search radius of parallax $\pi>3$ mas, limiting our sample
to stars at distances comparable to Taurus and Sco-Cen, two of the nearest
young associations. Both Taurus and Sco-Cen are thought to extend to a maximum
distance of approximately 200 pc (Preibisch & Mamajek, 2008; Galli et al.,
2019), so a minimum parallax of 3 mas (d$<$333 pc) allows for the exploration
of adjoining structures that may exist up to 100 pc beyond their known
extents, in addition to numerous other structures, both known and unknown.
This radius also includes nearly all of Perseus OB2 and the near edge of Orion
(Zari et al., 2019; Bally et al., 2008), allowing for partial coverage in
these more distant young associations.

We then impose several restrictions based on quality indicators, following
Arenou et al. (2018). Firstly we restrict the sample using the Unit Weight
Error ($u$), which can be interpreted as a goodness of fit measurement to the
astrometric solution in Gaia DR2. This value is given by
$u=\sqrt{\chi^{2}/\nu}$, where $\chi^{2}$ represents the $\chi^{2}$ value
between the source and a single-star astrometric
solution111astrometric_chi2_al in the gaia archive, and $\nu$ is defined as
the number of good CCD observations used in the astrometric
solution222astrometric_n_good_obs_al in the Gaia archive minus five (Lindegren
et al., 2018). We subsequently impose the following restriction on the Unit
Weight Error:

$\begin{multlined}u<1.2*\max[1,\exp(-0.2(G-19.5))]\end{multlined}u<1.2*\max[1,\exp(-0.2(G-19.5))]$
(1)

as suggested in Arenou et al. (2018).

In addition, we impose a restriction on the BP/RP flux excess
factor333phot_bp_rp_excess_factor in the Gaia archive $E$, which is defined as
the sum of the fluxes in the Gaia BP and RP photometric bands divided by the
flux in the G band. The configuration of the Gaia photometric bands implies
that $E$ should be slightly larger than 1 for stars with good photometric
measurements, however the BP and RP bands are vulnerable to external
contamination since they are based on integrated flux in a small field around
the star rather than the profile fitting used for the G band (Evans et al.,
2018). This contamination can manifest in the value for $E$, and contaminated
sources can therefore be removed using the following restriction:

$\begin{multlined}1.0+0.015(G_{BP}-G_{RP})^{2}<E<\\\
1.3+0.037(G_{BP}-G_{RP})^{2}\end{multlined}1.0+0.015(G_{BP}-G_{RP})^{2}<E<\\\
1.3+0.037(G_{BP}-G_{RP})^{2}$ (2)

This nearly matches with the restriction proposed in Arenou et al. (2018),
with a slight modification to the factor in the upper bound, where Arenou et
al. (2018) uses 0.06 rather than 0.037. We found that for known members of
Upper Sco and Taurus (Luhman, 2018; Rizzuto et al., 2015; Preibisch et al.,
2002, 1998), the slightly less restrictive conditions would occasionally show
very old photometric ages for these young stars. The more restrictive
condition removed nearly all such sources, significantly reducing our rates of
false negatives among known young objects.

We also impose restrictions on the number of visibility periods used in the
astrometric solution, which are defined as groups of observations separated by
at least four days. This requires that the astrometric solution is based on a
strong baseline of measurements.

$\text{visibility\\_periods\\_used}>8$ (3)

Finally, we restrict the sample based on the parallax inverse fractional
error444parallax_over_error in the Gaia archive ($\pi/\sigma_{\pi}$) to
exclude stars with a poorly constrained distance

$\pi/\sigma_{\pi}>10$ (4)

While the combination of cuts employed here excludes nearly half of known
young stars in Taurus and Upper Sco, these restrictions are important to
ensure that all absolute photometric measurements can be considered
trustworthy, which is critical to the reliable derivation of stellar
properties, particularly age. Once stars that pass these rigorous quality cuts
are used to identify and define the extents of young groups and associations
in space-velocity coordinates, the restrictions can be relaxed in future
targeted studies, which can significantly increase the completeness. See
Section 3.2 for a more complete and quantitative analysis into the
completeness of our results.

##  3 Methods

The solar neighborhood contains vast quantities of stars, including nearly 5
million within our 333 pc search limit that survive our Gaia quality cuts.
However, assuming a constant star formation rate over the last $\sim$10 Gyr,
only about 1 in 200 stars will have formed in the last 50 Myr. The Sco-Cen
association is thought to dominate nearby young populations, and it contains a
membership of only $\sim$10000 stars covering a massive 80 degree-long swath
along the plane of the Milky Way (Preibisch et al., 2002; Rizzuto et al.,
2015), so field populations will dominate young structures in most cases. To
isolate young populations from the field, photometry can often be used, as
many low-mass young stars have yet to leave the Hayashi track (Hayashi, 1961)
or have not yet settled onto the main sequence. Similarly, O and B stars are
readily identified as young, as they do not live long enough to become old.
However, in order to reliably identify young stars by photometric means, it is
critical that we understand what young stellar populations can and cannot be
confidently separated from the contamination of older field stars.

In this section, we describe our approach for identifying and characterizing
young stars in the Solar neighborhood. In Section 3.1, we describe the
implementation of a Bayesian statistical classification approach that uses a
model stellar population to identify credibly young stars in the solar
neighborhood, and we assess the success of this method in Section 3.2. The
rest of this section describes the methods used to further characterize stars
and the larger groups they might belong to. Section 3.3 describes our method
to obtain age estimates for individual stars, while we identify groups and
other young structures in the Solar neighborhood using the HDBSCAN clustering
algorithm in Section 3.4, lastly computing more precise bulk ages for these
larger groups in Section 3.5.

###  3.1 Generating the Model Population

With the excellent photometric and astrometric data provided by Gaia, we are
able to precisely determine the absolute magnitudes and colors of stars in the
sample. However, complicating factors such as metallicity, multiplicity, and
reddening modify the location of a star in absolute magnitude space, meaning
that a star that appears to be photometrically young may be better explained
by some combination of other factors (e.g., Sullivan & Kraus, submitted). We
therefore need to generate posterior distributions in age for each star and
marginalize over the other factors that may influence the photometric youth of
a star, including reddening and extinction, multiplicity, metallicity, and
stellar mass. To do this, we generate a simulated population of Each . We then
generate posterior distributions in age, as well as mass, corresponding to the
location of real stars from Gaia DR2 in the model population. By integrating
the age posterior over all ages less than 50 Myr, we can estimate the
probability $P_{Age<50Myr}$ that a star is genuinely young, allowing us to
isolate young stellar populations in the solar neighborhood.

####  3.1.1 Age

For the purposes of developing an age distribution, we can approximate the
star formation in the solar neighborhood to be constant over its lifetime, for
which we take the age of 11.2 Gyr from Binney et al. (2000). This is a
simplification of the true star formation history of the solar neighborhood,
however, as multiple studies have identified bursts in star formation activity
on scales of $\sim$ 1 Gyr (e.g. Rowell, 2013; Isern, 2019). Current literature
also indicates the possible presence of larger-scale trends in the star
formation rate, with most asserting some form of a decrease in star formation
since the galaxy’s earliest star formation bursts (e.g. Schönrich & Binney,
2009; Aumer & Binney, 2009). Most studies that assert a decrease in the star
formation rate take into account dynamical heating, a process by which older
stars are raised higher in the disk, therefore excluding many of them from
consideration. This effect works in opposition to a decreasing star formation
rate, implying that despite a possible gradual drop in the star formation
rate, the stars that remain within our Solar neighborhood are closer to having
an even distribution in ages.

We also assume that the age distribution is spatially isotropic, which allows
us to generate a single representative population that can be applied to all
of our Gaia stars, regardless of location. However, most recent star formation
is thought to have occurred within $\sim$100-200 pc of the galactic plane
(Urquhart et al., 2014; Anderson et al., 2019), with stars higher above the
plane mainly consisting mainly of older objects, potentially raised into their
current location by dynamical heating. As such, young stars will represent a
significantly higher fraction of the total stars close to the disk rather than
at high galactic latitudes. The assumption of an isotropic age distribution
may therefore inflate the rate of false positive young star identifications
near the Galactic poles by virtue of the absence of real young stars there,
but our rate of recovery among genuinely young stars should be roughly
constant across the entire population.

Not all literature agrees with the conclusion of a SFR that is either
decreasing or steady on average (e.g. Rocha-Pinto et al., 2000), however the
deviations that these alternative star formation histories might impose would
be of order unity and therefore have a limited impact on our posterior
distributions. On the main sequence, the photometry of stars is relatively
constant, so there is usually little photometric difference between otherwise
identical stars at, for example, 1 Gyr and 10 Gyr. Our objective in generating
the population of older stars should therefore be to ensure that their numbers
are approximately in proportion with their abundance in the solar
neighborhood, meaning that any star formation bursts can be essentially
averaged over. Any larger-scale SFR trends may slightly affect our results by
means of changing the relative fraction of young and old stars, however these
effects can be accounted for by modifying the probability cutoff for selecting
candidate young stars in Section 3.1.7.

A star formation rate that is to first order uniform over the last $\sim$11
Gyr implies that stars formed within the last 50 Myr should account for less
than one percent of the total stars in the solar neighborhood. Subsequently, a
representative sample of ages would overwhelmingly skew towards very old
objects, leaving a sparse sample of young stars. Since stars evolve quickly on
the pre-main sequence, a sparse sample of model stars at these young ages
would be insufficient to effectively assess whether any observed star is
genuinely young. We therefore draw model stellar ages from a Log10-uniform
distribution spanning 1 Myr to 11.2 Gyr. This does introduce an exponential
bias towards young objects, however that can be corrected using a prior, as
described in Section 3.1.7.

####  3.1.2 Primary Masses

The initial mass function, which displays how the abundance of stars varies
with mass, is the basis for this project’s mass generation. It follows a
distribution that has been gradually revised since the first form presented by
Salpeter (1955), with the most recent revisions from Chabrier (2003) and
Chabrier (2005) using a log-normal distribution at low mass (M $<$ 1 M⊙), and
the Salpeter (1955) power law distribution for higher masses. The IMF has been
consistently shown to be invariant based on environment, so spatial variations
to our mass generation are not of concern (Chabrier, 2005; Offner et al.,
2014). Slightly different log-normal solutions have been presented for
individual stars and complete stellar systems, and since we require masses of
primaries that may or may not be in multiple systems, we make use of the
Chabrier (2005) individual star Initial Mass Function (IMF) to generate the
masses of primaries. Both the individual and system Chabrier (2005) initial
mass functions are plotted in Figure 1, alongside the distribution of system
masses we reach after adding stellar companions in Section 3.1.4.

The minimum mass we consider is 0.11 M⊙, which is roughly the smallest mass
that is available in all PARSEC v1.2S isochrones (Chen et al., 2015)). Young
stars of this mass also reach Gaia’s approximate limiting magnitude of
G$\simeq$21 near our maximum distance of 333 pc (Arenou et al., 2018), so
objects in this mass range will be accessible throughout our entire search
radius. Our maximum mass was set to 20 M⊙, which excludes a small subset of
stars that are both very rare and very bright. Main sequence stars of more
than 20 solar masses have absolute magnitude G$<$-3 (Chen et al., 2015), which
corresponds to a magnitude of roughly 3 at the distance of Sco-Cen. Gaia
photometry and astrometry in this range is known to be poor due to the start
of pixel saturation at G$\sim$6 (Arenou et al., 2018), so the contribution
Gaia can make towards furthering the understanding of these very bright stars
is relatively minimal.

Due to a combination of the higher average mass of primaries relative to the
complete population considered in the individual IMF, as well as the mass
limit of 0.11 M⊙ imposed on all model stars, our system mass distribution does
not quite replicate the Chabrier (2005) system IMF. For stellar systems, which
are discussed in depth in Section 3.1.4, a minimum component mass of 0.11 M⊙
implies a minimum binary mass is 0.22 M⊙, and as a result this mass marks a
small local minimum in the mass function in Figure 1. As such, the
discrepancies that do exist are mostly brought on by limitations in the masses
of objects our system can create. However, despite these minor deviations, the
final simulated system mass function remains comparable to the system IMF also
presented in Chabrier (2005), never differing by more than 30%. We expect the
effects of a slightly discrepant model IMF to be relatively minor in our final
results. The overabundance of higher mass stars relative to lower mass
binaries may overweight the contamination from subgiants, however,
particularly in the parameter space occupied by the pre-main sequence, G and F
stars follow evolutionary tracks well-separated from the early M stars, which
show the greatest underabundances in our model (Chen et al., 2015). This
suggests that situations where the relative abundances of these different mass
model stars has a significant impact on our results will be very rare. Any
such effects will be further minimized by our generous 10% probability
threshhold for the consideration of a star as young, and the reliance on
kinematics to cull likely erroneously-detected field stars, as outlined in
Section 3.4.

Figure 1: The Chabrier (2005) individual object Initial Mass Function (bright
yellow), which was used to set the primary masses in our model systems,
plotted against the distribution of system masses in our model after binary
and triple companions are added (green). The Chabrier (2005) system IMF is
included for reference (dark yellow). Our model systems never deviate from the
Chabrier (2005) system IMF by more than 30%, although our model does show a
minor deficiency of low-mass binaries with components below the lower mass
limit of our model grid (Chen et al., 2015).

####  3.1.3 Metallicity

While stars are thought to be more metal-rich on average at later formation
epochs, the age-metallicity relationship in the solar neighborhood appears to
be relatively flat beyond an initial early enrichment phase (e.g., Lin et al.,
2020; Haywood et al., 2019), inferring that a single non-age-dependent
metallicity distribution will be sufficient to represent the solar
neighborhood. Regardless of whether metallicity is truly constant with time,
it may not influence our detection of young stars, since the nearly age-
independent photometry of main-sequence stars means that the age will weigh
minimally on the photometry of these stars, regardless of metallicity.

A possibly more consequential bias that a single metallicity distribution may
introduce relates specifically to young stars, as most that are known have
solar metallicities. Using the full range of possible metallicities to
represent these stars may therefore present a wider range of possible
photometry compared to what actually exists, particularly for low
metallicities. This would result in an underestimation of $P_{Age<50Myr}$, as
these low-metallicity stars have photometry more consistent with older
populations (Chen et al., 2015). However, since we wish to identify structures
beyond those already known in the solar neighborhood, we have little cause to
tether our populations to an assumption of solar metallicity. Even if the
assumption of near-solar metallicity is universally appropriate, placing stars
too low on the pre-main sequence due to low metallicity is a conservative
error, and will therefore not lead to spurious assessments of youth. The
resulting lower probability estimations can subsequently be negated using more
permissive probability cuts for the identification of young stars. We
therefore conclude that a non-age dependent distribution for metallicity would
provide a reasonable approximation to reality. Data Release 2 of the GALAH
survey (Buder et al., 2018; Hayden et al., 2019) provides an empirical sample
containing over 62000 stars within 500 pc of the sun which we use as a
representative sample for the metallicities of stars in the solar
neighborhood.

Our metallicities were randomly generated from a slightly modified version of
the GALAH metallicity distribution. We first restricted the metallicities to
between [Fe/H] = -1.0 and 0.5. The lower limit was due to the negligible
number of GALAH stars with metallicities in that range ($\sim$ 0.1%), while
the upper limit was chosen to comply with the metallicity range available in
PARSEC v1.2S isochrones: -2.2$<$[Fe/H]$<$0.5 (Chen et al., 2015). There are
more stars that exceed this high metallicity limit in the GALAH data compared
to the low-metallicity limit, however the fraction of total stars occupied by
these outliers remains under 0.5% (Hayden et al., 2019). Therefore, these
extremely metal-rich stars are not sufficiently numerous for their exclusion
to dramatically skew our results. We then binned the GALAH measurements with
[Fe/H] between -1.0 and 0.5 into 200 evenly-spaced bins over that range, and
smoothed the resulting distribution using a Savitsky-Golay filter with a
15-bin window length and third-order polynomial. This generated a smooth
metallicity distribution consistent with the original histogram from which we
can draw samples, which is shown in Figure 2

Figure 2: The metallicity histogram for stars observed by the GALAH survey.
The dark curve shows the smoothed distribution that we generated metallicities
from (prior to normalization).

####  3.1.4 Multiplicity

Binaries are quite common in the Milky Way, with a mass-dependent abundance
fraction that increases in tandem with the mass of the primary (e.g., Duchêne
& Kraus, 2013). For stars less massive than the sun, between one fifth and one
half of all stars are binaries, and depending on the distance many of these
cannot be resolved by Gaia. While resolved binaries are not of concern for the
purposes of this study, as Gaia is able to generate accurate photometry and
astrometry for each resolved component, unresolved binaries create combined
photometric detections that appear overluminous relative to comparable single
stars. Dim companions have a relatively minor impact on the system photometry,
however the presence of an unresolved equal-brightness binary can double the
incident flux and subsequently lower the magnitude by up to 0.76 mag, a much
more significant contribution. This brightness relative to the main sequence
is comparable to what we see on the pre-main sequence, and therefore these
older unresolved binaries can easily be mistaken for younger single stars.
Therefore, to accurately assess whether a star is young, the abundance of
multiple stars, distribution of companion masses, and companion separations
must all be well-modelled to quantify the probability that stars in any region
above the pre-main sequence are multiple systems rather than young stars.

We therefore add multiples to our synthetic population, considering only
single, double, and triple systems to avoid complications with the
architecture of rare, higher-order systems. We begin with the statistics
reported in Duchêne & Kraus (2013), including the multiplicity frequencies
(MF), which measure the occurrence frequency of multiple systems of any kind,
and companion frequencies (CF), which considers the number of companions per
primary star (originally sourced from Raghavan et al., 2010; Delfosse et al.,
2004; Dieterich et al., 2012; Kouwenhoven et al., 2005, 2007; Abt et al.,
1990; Sana et al., 2012; Mason et al., 2009; Chini et al., 2012; Preibisch et
al., 1999). While the MF describes the probability of a star being in some
form of multiple system, the CF allows for conclusions to be drawn about the
occurrence frequencies of both double and triple stars.

To generate a smooth distribution, we fit a logarithmic function as in
equations 5 and 6 to the MF and CF values from Duchêne & Kraus (2013), using
reasonable intermediate values for the MF in cases where only upper limits are
provided. These fits are shown in Figure 3, matching the linear trends in
semi-log mass vs MF and CF that are followed by the Duchêne & Kraus (2013)
values for MF and CF. The resulting solutions for CF and MF as a function of
primary mass are as follows:

$MF=0.272*\text{log}_{10}M+0.431$ (5)

$CF=0.491*\text{log}_{10}M+0.625$ (6)

Then, assuming that the probability of each multiplicity follows some
geometric series $Ck^{m}$ (as in Duchêne & Kraus, 2013), where $m$ is the
multiplicity of the stellar system, we fit for the unique solution for C and k
that returns the appropriate values for CF and MF interpolated from our fit.
The resulting abundance fraction distributions for multiples, binaries, and
higher order systems are given in Figure 4. The result has the correct edge
behavior, as the double and especially triple rates drop dramatically as the
primary mass approaches our minimum mass limit, and the multiple rate
approaches but does not exceed one at the high-mass limit.

The probability distribution of the resulting companion masses is typically
described as a function of the mass ratio between the primary and secondary
components $m_{s}/m_{p}=q$, following a power law distribution of the form
$P(q)=Cq^{\gamma}$ (Kraus et al., 2011; Rizzuto et al., 2013). We therefore
generate secondary and tertiary components based on the mass ratio probability
distributions presented in Kraus et al. (2011) and Rizzuto et al. (2013),
which cover $q$ values between 0.1 and 1. In order to exclude companions in
the brown dwarf regime that are not massive enough to be included in the
PARSEC isochrones, we multiply these probabilities by an additional factor,
which limits the section of the mass ratio distribution from which companions
can be drawn to only that region between the primary mass and the minimum mass
accessible through the PARSEC isochrones. The corrective factor for a given
primary mass is therefore the integral of the probability distribution over
the window from the primary mass to the minimum PARSEC mass, normalized such
that the sum over all $q$ is 1. Once this corrective factor to the
multiplicity probabilities is applied, the primaries are separated into three
populations, each with distributions for $q$ governed by a power law with a
different exponent $\gamma$: $\gamma=0.4$ for M $<$ 0.7 M⊙, $\gamma=0$ for
0.7-2.5 M⊙, and $\gamma=-0.46$ for masses exceeding 2.5 M⊙. The two lower mass
ranges come from Kraus et al. (2011), and the upper mass bracket comes from
Rizzuto et al. (2013).

The generation of masses for all primaries and companions represents the final
step necessary before computing simulated inherent photometric properties for
these stars. The mass distribution of systems after the inclusion of
companions is shown in Fig. 1, plotted against the Chabrier (2005) IMFs.

Figure 3: The fits to the MF (blue) and CF (red) distributions as a function
of primary mass. Note that the confidence intervals for the two rightmost
points for the MF distribution are lower limits, and the dot provided is the
chosen intermediate value. These points are heavily de-weighted compared to
the well-constrained MF values at lower masses, and have minimal impact since
our selected sample consists almost entirely of stars of spectral type G or
later. Figure 4: The final distributions used for the binary, higher order
(which are treated as triples), and total multiplicity fractions. The blue
line is equivalent to the MF fit in Figure 3.

The separation of multiple system components must also be considered, as only
spatially unresolved systems will display merged photometry. The stellar
populations we wish to draw young stars from come from the Gaia mission, and
since Gaia reports resolved stars independently, multiples that are resolved
in Gaia are seen as a series of two or more singles. To generate approximate
separations, we take the distribution of periods presented in Raghavan et al.
(2010), randomly select an orbital period from that distribution, and compute
an orbital semi-major axis from the mass of the system using Kepler’s Third
Law. A corrective factor of $\pi/4$ is applied to account for the projection
of these orbits on the plane of the sky, which is modulated by the presence of
randomly-oriented systems. Recent literature investigating the separations of
M dwarf binaries has shown tighter lognormal distributions centered at
slightly smaller periods for these stars (e.g., Bergfors et al., 2010),
however these results are not currently well-enough constrained for us to
justify including a mass-specific period distribution. We may introduce
revised separation distributions in a future iteration of this work, as the
current implementation likely underestimates the abundance of unresolved
binaries for M dwarfs. However, even if an updated form of the separation
distribution were to imply that all M-dwarf binaries are unresolved, this
would only represent a roughly 30% increase to the unresolved rate for
companions to a $\sim$0.2 M⊙ star, so we do not expect significant effects
from the use of the Raghavan et al. (2010) period distribution.

Most Gaia stars with separations exceeding 1 arcsecond are resolved as
independent sources (Rizzuto et al., 2018; Ziegler et al., 2018), although
this does not appear to be a hard boundary, as there is likely a population of
systems at separations below 1 arcsecond that are neither unresolved nor
reported independently, a factor which may reduce the mean separation at which
this transition to independent identification occurs. However, even with that
transition separation lowered to 0.2 arcseconds, the unresolved rate for an
average 1 M⊙ star only drops by 25%, inferring that the use of 1 arcsecond as
a hard limit for the transition between unresolved and resolved populations
appears to be a reasonable approximation. We take a distance of 234 pc, which
is the mean distance to the Gaia stars investigated, and break up test stars
with separations exceeding 1 arcsecond at that distance. Stars with smaller
separations are marked as unresolved systems, which are combined
photometrically into a single star. For triple stars, separations are
generated for each companion. If both are unresolved, then all three members
are photometrically merged. If only one model separation is resolved, then the
second component is counted as resolved from the primary, and the third is
randomly assigned to either the primary or secondary. The result is that the
system is now treated in our model as one unresolved double with a resolved
single nearby. The choice to separate binaries at a single fixed distance,
rather than make it adaptive to the distance of the target star, was an
approximation made to reduce computational costs, however we do not expect the
effects of this choice to be significant. At 10 pc, a companion to a 1 M⊙ star
would have an approximately 34% chance of being unresolved, compared to about
72% for a star at our field limit of 333 pc, and 68% for our chosen distance
at 234 pc. As such, the only major inaccuracies will be at very small
distances, where we already expect our cluster detection to be less sensitive.
Within the nearest 50-100 pc, geometric effects on the 2-d Gaia transverse
velocities make it increasingly difficult to detect stars with common 3-d
motions, making structures difficult to identify regardless of the recovery
rate of their stars (see Section 3.4).

####  3.1.5 Generating Photometry of Simulated Stars and Systems

For each mass, metallicity, and age generated, we use isochrones to create
corresponding Gaia G, GRP, and GBP magnitudes. We make use of PARSEC version
1.2S555http://stev.oapd.inaf.it/cgi-bin/cmd isochrones made freely available
online (see Tang et al., 2014; Chen et al., 2015; Bressan et al., 2012).

Note that these isochrones do not include the white dwarf cooling sequence,
and therefore these objects are treated as “dead”, with no photometric
contribution. Since the photometry of these objects says little about the
white dwarf’s age relative to the formation of the system, this exclusion
should have a minimal impact on the results. These objects also occupy a very
different parameter space on the CMD, so they are easy to remove and nearly
impossible to mistake for young main sequence stars.

The isochrones cover a metallicity range evenly spaced between [Fe/H] = -1.0
to 0.5, and a less regular age selection, which was designed to optimize the
sampling of parameter space between the ages from 6$<$log10(age/yr)$<$10.049.
The main concern when generating a grid to draw our synthetic population from
is that with a sufficiently sparse sample in age, stars of a given mass would
not have samples in between, for example, the main sequence and giant branch.
This can result in situations where any interpolation on these photometric
grids can end up deriving photometry somewhere between the main sequence and
giant branch, which may not follow the evolutionary track stars follow between
those sequences. To minimize artefacts caused by under-sampling, we require
that for each slice of the grid in both mass and age, at least two points in
age populate the horizontal branch, except for the most extreme massive stars.
This requirement ensures that at least the lower red giant branch is well-
covered by our grid, with the top of the RGB and Horizontal Branch being
somewhat less well-covered. Further grid resolution beyond covering the lower
RGB would only increase the computational cost without improving our
calculations, as these more evolved stars do not overlap with the pre-main
sequence stars in photometric parameter space.

The PARSEC v1.2S isochrones give a list of masses for each metallicity/age
combination, and we regrid these onto a grid with 0.11$<$log(M/M⊙)$<$20) and
non-uniform point density. For masses, very sparse sampling is sufficient for
stars less massive than the main sequence turnoff at the end of our age
bracket, 11.2 Gyr. All stars below this point do not ever evolve off the main
sequence, and therefore the only quick evolutionary process with sensitivity
to mass that must be captured is the descent along the Hayashi track, which is
significantly slower than evolution along the RGB and much less mass-sensitive
(Hayashi, 1961; Chen et al., 2015). Not many grid points are required to
capture that evolution in the mass dimension, so only 100 linearly-spaced grid
points are used for masses M $\lesssim$ 0.8 M⊙. We then use 4000 bins between
0.8 M⊙ and 2.5 M⊙, another 4000 between 2.5 M⊙ and 10 M⊙, and 600 more bins
between 10 M⊙ and 20 M⊙, for a total of 8700 bins. These binning choices do
result in the under-sampling of the more evolved stars in the population,
however the pre-main sequence and stars sharing that parameter space are all
well-covered.

We assigned model stars random ages, metallicities, and masses according to
the relevant distributions, and linearly interpolated magnitudes for each star
from our isochrone grid. When identified as part of an unresolved binary or
triple system, we added the magnitudes of the component stars together,
creating a merged photometric result for each unresolved system. This produced
intrinsic synthetic G, BP, and RP magnitudes for any model star or system,
resulting in the color-magnitude diagram presented in Figure 5. While this
completes the generation of our model population, these intrinsic magnitudes
have yet to be modulated by the conditions in which a real Gaia star might
exist, most notably extinction, which we address in the next subsection.

Figure 5: A color-magnitude diagram for the stars generated by simulating a
population with our randomly generated ages, masses, metallicities, and
stellar multiplicities, then computing observables using the PARSEC
isochrones.

####  3.1.6 Extinction

Significant interstellar extinction moves the main sequence of a stellar
population on the HR-diagram to the right and downwards of where it would be
without that extinction. Consequently, heavily extincted stars can appear
younger than they would without the effects of extinction. Due to the
significant regional variations in extinction, the use of a single,
directionally independent distribution to modulate our simulated photometry
would fail to differentiate between young stars and those that are merely
heavily extincted. Instead, we make use of three dimensional, all-sky
extinction maps, and generate a unique reddening distribution for each Gaia
star based on that map.

We use the STILISM maps of reddening in the solar neighborhood from Lallement
et al. (2019), which use Gaia DR2 and 2MASS photometry and distances to
compute reddening for a large volume centered on the solar system. We
interpolate reddening parameters for each Gaia star based on its on-sky
location using a Monte-Carlo framework with distances drawn from the Gaia
parallax and parallax uncertainty. This framework captures how Gaia distance
uncertainties impact reddening uncertainties.

An additional uncertainty modulation is added to each interpolated reddening
result to reflect the inherent uncertainty in the maps themselves. Lallement
et al. (2019) reddenings often have very asymmetric uncertainties, so for each
generated reddening and corresponding uncertainties, we select a random value
from a normal distribution, multiplied by the uncertainty on that side of the
reddening value recorded in the STILISM maps. These uncertainty results are
added to each randomized reddening, allowing us to generate probability
histograms for the reddening of all stars in the Gaia set. The reddenings of
the simulated stars are therefore drawn from these histograms, separately for
each Gaia star. The relationships presented in Wang & Chen (2019) were used to
relate E(B-V) to extinction in the Gaia BP, RP, G filters. While these
conversions are expected to differ somewhat by the effective temperature of
the target, these effects will be negligible compared to the intrinsic sources
of uncertainty in the Lallement et al. (2019) reddening maps, especially given
that the Wang & Chen (2019) reddening conversions are based on relatively cool
red clump stars which are photometrically similar to many of the pre-main
sequence stars we focus on in this paper. The addition of reddening completes
a population of simulated stars with a reddening distribution reflecting the
reddening expected of each Gaia star based on its location.

####  3.1.7 Generating Star Statistics

The completed population of sample stars allows us to generate statistics for
our population of Gaia DR2 stars. For each real star, we use Bayes’ theorem to
generate the probability that each sample star is consistent with the Gaia
photometry according to the following formulation, as implemented in Huber et
al. (2016):

$p(y|g)\propto
p(y)p(g|y)=\prod_{i}exp\left(-\frac{(g_{i}-x_{i})^{2}}{2\sigma_{g,i}^{2}}\right)$
(7)

where y is a sample star with inherent properties $y$ = {age, [Fe/H], mass,
E(BP-RP), multiplicity} and with observables $x$ = {G, BP, RP}. We then
compare the observables of our sample stars to our Gaia star $g$ with the same
observables, converting from apparent Gaia magnitude to absolute magnitude
using the Gaia parallax measurements.

These distances were generated by inverting the parallaxes, and while this is
an imperfect estimate for high-relative uncertainty parallaxes, the distance
estimates from Bailer-Jones et al. (2018) follow the inverted parallax very
closely for the stars with $\pi/\sigma_{\pi}>10$ that we investigate here. The
subscript $i$ multiplies over each observable in x, and $\sigma_{g,i}$ is the
uncertainty in the Gaia measurement of the observable $i$, including both
photometric and parallax uncertainties.

In most cases, the prior $p(y)=1$, owing to the consistency of our
metallicity, mass, reddening, and binarity rates with true values. Age,
however, does require a non-trivial prior, as we significantly over-sample
young stars. We therefore add a prior equal to the age of the star, which
nullifies the exponential bias cause by the uniform sampling in log-space,
producing a flat probability distribution in age for the sample stars. To
produce any probability distribution, we simply sum the $p(y|g)$ values for
each model star into a histogram, and normalize over all sample stars. The
probability of any given condition on a property, such as $P_{Age<50Myr}$ can
therefore be found by summing over the relevant bins. The resulting
$P_{Age<50Myr}$ values for each star in the Gaia sample are presented through
the color-magnitude diagram in Figure 6. Most stars have $P_{Age<50Myr}$
similar to their neighbors, with some exceptions in cases with larger
uncertainties, or cases where either the BP or RP filter is discrepant, which
may not translate well to this 2-d plot from the 3-d input magnitude data set.

All Gaia stars with less than 50 simulated stars within 1-sigma of them in
color-magnitude space were culled from the sample, as we identify them as
inconsistent with the stars in our simulated population. As we do not consider
white dwarfs and white dwarf-main sequence binaries in our model, most of
these stars are excluded, as is evident from the near-absence of a white dwarf
cooling sequence in Figure 6. To extract a population of young stars from our
sample, all objects with $P_{Age<50Myr}<0.1$ were removed, leaving a
population of 30518 credible young star candidates out of the full sample of
nearly 5 million Gaia sources. These candidate young stars are presented in XY
galactic spatial coordinates in Figure 7, and compiled into a master candidate
list in Table 1. As a result of the reddening vector direction shown in Figure
6, the underestimation of reddening, especially near the subgiant branch, does
have the ability to make older stars appear young. Minor anomalous reddening-
related clumps do appear. While we do detect minor spatially clustered but
kinematically scattered clumps of stars that appear consistent with anomalies
from local reddening underestimates, many more young structures are visually
identifiable, such as Sco-Cen, Orion, Perseus, and Taurus.

Figure 6: A color-magnitude diagram containing all Gaia stars with more than
50 simulated stars within 1-sigma of their locations in the CMD. All
magnitudes are absolute and dereddened. The colors represent the probability
that the star is young, defined by this paper as younger than 50 Myr. The
reddening vector for AG=2 is labelled, showing the direction along which older
stars may move to interfere with the pre-main sequence if not dereddened.
Figure 7: All credible young stars identified by our search to have
$P_{Age<50Myr}>0.1$, plotted in X/Y galactic coordinates. The circled region
is an example of a reddening anomaly, where the underestimated reddening
causes stars to be marked as young. That particular region is in the direction
of heavy reddening in Aquila and Serpens. This group is not identified in the
clustering stage (see Fig. 10), as the stars do not have common motions, and
are likely dominated by reddened field stars.

###  3.2 Recovery

The method we employ in this paper searches the solar neighborhood for stars
with a significant probability (P$>$0.1) of an age less than 50 Myr, and we
should therefore have demonstrable capabilities to identify stars in that age
range. Here we investigate the efficiency with which those nearby young stars
are identified, and what kinds of stars are detectable for groups with
different ages. There are two main ways in which stars might be missed in this
survey: failure of Gaia photometric or astrometric quality cuts and
misidentification by our pipeline as old. In this section we quantify the
losses due to each of these factors, and provide insight into why some of
these losses take place.

The Gaia quality cuts, which remove stars based on factors like goodness of
fit to the Gaia photometric and astrometric models, parallax error, and the
number of Gaia visibility periods used (see Section 2.1 for a full explanation
of these) are necessary to ensure that our Gaia stars have reliable photometry
and astrometry. However, it is nonetheless useful to investigate the
properties of the stars lost to reveal any potential biases the cuts might
impose. By comparing the samples of stars that pass and fail these quality
cuts, we can reveal what properties are most likely to lead to the rejection
of a star based on quality. We use a large sample of Taurus members drawn from
(Krolikowski, submitted) to make this comparison, which has a wealth of
information on $\sim$500 stars in the region, including spectral type,
extinction, and binarity. All of these factors may influence the quality of
any Gaia astrometric or photometric solutions. We find that our Gaia quality
culling appears to correlate very little with stellar spectral type over the
late G to M range, however those cuts do correlate strongly with binarity and
extinction. Unresolved binaries are almost twice as likely to be removed by
Gaia quality restrictions, with 59% of singles passing these cuts, compared to
only 31% of multiples. Binaries are known to both introduce additional
astrometric noise and skew the photometric solution (e.g., Arenou et al.,
2018), so their more frequent failure of quality cuts is expected. For heavily
reddened populations with $A_{V}>3$, only 20% of sources pass the quality
cuts, compared to 67% for $A_{V}<0.5$. This heavy loss rate in highly reddened
locations is consistent with our results later in this paper, which include
visibly incomplete samples of young stars in active star-forming sites such as
Perseus and Lupus.

Figure 8: Recovery rates after the imposition of Gaia quality cuts in Taurus
(green squares), Upper Sco (purple diamonds), $\beta$ Pic (orange circles),
and Tucana-Horologium (red triangles), binned by spectral type. In all bins,
the order of the groups in recovery rate matches that of their relative ages,
with the essentially newborn Taurus association having the most complete
recovery. The error bars display the binomial 68% confidence interval.

Next we investigate recovery rates from our pipeline for the identification of
young stars, which include only objects that pass the Gaia quality cuts. To
calculate recovery rates, we gather lists of known members for four well-known
associations: Upper Sco (in Sco-Cen), Taurus, the $\beta$ Pictoris Moving
Group, and the Tucana-Horologium Moving Group. These groups cover nearly the
full range of detectable cluster ages, spanning from Taurus, which is nearly
newborn ($<$ 5 Myr; Kraus & Hillenbrand, 2009), to Tucana-Horologium, which
has an age of approximately 45 Myr (Kraus et al., 2014; Bell et al., 2015). We
can therefore use these stars as a strong representative sample of the
populations we expect to identify in this paper.

As the youngest population we investigate in this section, stars in the Taurus
Association should be the easiest to identify reliably as young through purely
photometric methods. To verify this, we drew known members from the Esplin &
Luhman (2019) Taurus catalog, which compiled verified members from literature
and introduced several new members. These members are focused on the young
central subgroups in Taurus (i.e. Greater Taurus groups 8-11 in Section 4.1),
making them a more homogeneous sample for comparison relative to the more
diverse sample from Krolikowski (submitted), which includes both the objects
in Esplin & Luhman (2019) and those in more peripheral regions, which may be
older. Esplin & Luhman (2019) identified members as candidates using proper
motions and photometry, and then confirmed them using spectroscopic
observations including Li I absorption. We found that 189 of 351 (54%) stars
in Taurus have Gaia detections that survive the quality cuts described in
Section 2.1. While these restrictions did result in the loss of nearly half of
the members, it is important for the reliability of our results that objects
with poor Gaia data are excluded. Of the stars that survived the quality cuts,
172 of 189 were identified as young, or 91%. We compile the recovery rates as
a function of spectral class for Taurus and the other three test associations
in Fig. 8, where we show that the recovery rate in Taurus is essentially
complete for spectral classes M and K, with a sharp sensitivity decline for
G-type stars. This lowered sensitivity is caused by the photometric overlap
between pre-main sequence stars earlier than about G8 and the subgiant branch
on the HR diagram, which significantly reduces the probability that a star
there is young.

The Upper Sco Association is slightly older than Taurus (age $\sim$5-11 Myr;
Rizzuto et al. 2015; Pecaut et al. 2012; Preibisch et al. 1999), and we
therefore expect less sensitivity to early type stars as members begin to
merge onto the main sequence. The Upper Sco sample comes from Rizzuto et al.
(2015), Preibisch et al. (2002), and Preibisch et al. (1998), making use of
spectroscopic observations of Li and H$\alpha$ combined with complete
kinematic solutions to confirm the membership of the stars (Rizzuto et al.,
2015). A total of 318 members out of the total list of 477, or 67%, had a
high-quality Gaia counterpart that passed our quality cuts. Of those high-
quality Gaia detections of known Upper Sco members, 274 were identified as
young by our method, or 86%. As presented in Figure 8, our completeness is
best for the lowest-mass stars, which provide essentially complete results for
M dwarfs up to about M1. For stars earlier than M1, the sensitivity drops
gradually, with the same significant sensitivity drop we observed in Taurus
appearing beyond spectral type G8.

The older $\beta$ Pictoris Moving Group (age $\sim$ 23 Myr; Mamajek & Bell
2014), hosts a more elusive population, as these stars sit on a pre-main
sequence that is less well removed from the photometry of background stars. To
check the recovery rate in $\beta$ Pic, we used the catalog from Shkolnik et
al. (2018), which identified members using both 3-d motions and the presence
of Li and H$\alpha$. 94 of 173 stars in the association, or 54%, survived our
Gaia quality cuts. The rate at which we recovered these young stars was
considerably lower, however, with 44 out of 94 members being recovered, or
47%. Rather than there being a sensitivity limit around G8, the brightest
stars identified in $\beta$ Pic have spectral types in the early K, caused
primarily by the more luminous stars having already settled close to the main
sequence. Later K stars showed low recovery rates, suggesting that they are
usually only marginally identified (see Fig 8). Significantly better recovery
rates were observed for mid-M stars, peaking at approximately 80% for M3-M4,
suggesting that our method remains effective at consistently identifying the
less massive stars that remain well above the main sequence even at older
ages.

The Tucana-Horologium association is near the 50 Myr target age limit with an
age of $\sim$45 Myr (Kraus et al., 2014), and therefore we should expect
significantly lower sensitivity to this population. The sample population of
Tuc-Hor members that we used comes from Kraus et al. (2014), which, like for
our $\beta$ Pic population, used a combination of 3-d kinematic information
and stellar youth indicators to verify membership. 80% of Tuc-Hor members (or
115 out of 143) survive the Gaia quality cuts. As expected, our recovery of
Tuc-Hor members showed a significant decay in the recovery rate of stars
overall, with no noteworthy sensitivity for stars earlier than a spectral
class of about M1 (see Fig. 8). The recovery rate for the later M-dwarf bins
in Tuc-Hor is around 20%, which, while significantly weaker than our recovery
for younger stars, does still identify 17 Tuc-Hor members out of the complete
sample of 115 across all spectral types, or 15%, enough to make the population
potentially identifiable.

The recovery rates among stars with quality Gaia detections in each of the
four regions discussed here are compiled in Figure 8, showing clearly how our
recovery rates vary with spectral type, and with the age of the stars being
observed. There we show how the effectiveness of our method at a given
spectral class negatively correlates with the age of the region, with the
nearly newborn Taurus Association being essentially complete from G8 to late M
excluding the quality cuts, and the 45 Myr old Tuc-Hor Moving Group having
extremely limited sensitivity essentially restricted to M dwarfs. For ages in
between, sensitivity gradually falls off, with G and then K stars being lost
first, and M dwarf sensitivity beginning to drop past $\sim$ 20 Myr as for our
$\beta$ Pictoris sample. For Tuc-Hor, an association near our target age
limit, we maintain a recovery rate of approximately 20% among M stars, while
only one earlier star is identified. Despite the restrictions to our
sensitivity for these older groups, our ability to detect non-negligible
numbers of members in even the relatively old Tucana-Horologium Moving Group
suggests that large enough groups will remain identifiable up to the upper
edge of our target age range at 50 Myr.

Despite the successful stellar recovery among these older populations, the
composition of the stars recovered should be treated with caution, as we
expect stars with features that artificially inflate a star’s photometric
youth, such as unresolved binaries, to be detected much more easily than an
average star in the samples included in this section. With an overall recovery
rate around 15%, as for Tuc-Hor, unresolved binaries may begin to dominate the
recovered sample. Particularly in these older populations, demographic studies
can be greatly expanded by using the locus of an association in space-velocity
coordinates as defined by the limited sample we identify as young to locate
and reintroduce probable members that are not recognized as photometrically
young by our pipeline. The reintroduction of photometrically older probable
members can not only expand the completeness of the original sample, but also
help to suppress biases produced by the potential overrepresentation of
binaries, such as in the derivation of group ages, an application we explore
in Section 3.5.

###  3.3 Age Estimation

The primary objective of generating our model population and corresponding
posterior distributions was the identification of young stars in the presence
of diverse stellar populations in the solar neighborhood. However, the age
distributions that we generate can also be used to derive approximate ages for
most of the young stars in our sample. For candidate young stars, the age
probability distribution from our pipeline consists of a peak at a young age,
and a significant tail

Figure 9:

While binaries are very common, potentially representing a majority of young
systems (e.g. Kraus et al., 2011; Raghavan et al., 2010), their companions
introduce a very wide range of photometric excesses, depending on the mass of
the companion. As such, photometry is a much less consistent age diagnostic
for multiple systems compared to stars without a companion. We therefore
exclude binaries and higher-order systems from our individual age estimates,
subsequently assuming that all Gaia sources have no companion. For single
stars, the lack of consideration for binaries will ensure that our results are
not influenced by the wide range of photometric solutions possible from the
presence of an unresolved companion, making the age estimate of high quality.
If a star does have a companion, it will contribute to a background of stars
with younger age fits that are not consistent with one another, making them
possible to clip out when looking at populations.

We also fix the assumed metallicities of our candidate young stars to the
solar value. This assumption has been shown to be appropriate for the vast
majority of nearby young clusters and associations (e.g., Viana Almeida et
al., 2009; Mamajek, 2013), which have a uniformity in their metallicities
consistent with the thorough mixing of the ISM on 100pc-level scales (e.g. de
Avillez & Mac Low, 2002). Table 3 in Viana Almeida et al. (2009) presents
metallicity values for most nearby young associations, with all solutions
falling between [Fe/H]=-0.13 and [Fe/H]=+0.04, mostly within one standard
deviation of solar. The $\alpha$ Persei cluster is one of the most metal-rich
known groups accessible through our search, and it has a metallicity of only
[Fe/H]=+0.18 (Pöhnl & Paunzen, 2010). Due to the mixing of the local ISM, we
would not expect any of our groups to depart significantly from the trend of
consistently near-solar metallicities for young nearby associations and
clusters. Unlike these young groups, older clusters have a more diverse set of
metallicities reflective of those in the rest of the local stellar population
(Buder et al., 2018; Hayden et al., 2019, see Fig. 2). The Hyades and Praesepe
have notably super-solar metallicities at [Fe/H]=+0.146 and [Fe/H]=+0.21,
respectively (D’Orazi et al., 2020; Cummings et al., 2017, e.g.,), while
numerous other clusters are much more metal poor, such as M35 and NGC 2506 at
[Fe/H]=-0.21 and [Fe/H]=-0.52, respectively (Bouy et al., 2015; Friel & Janes,
1993). This contrast in the metallicity variation between the youngest stellar
populations and various older populations in the solar neighborhood
strengthens the decision to restrict the metallicity to solar for our in-depth
look at young stars, while assuming a much broader prior metallicity
distribution for our wider population.

To extract age estimates from our output age probability distributions, we
feed model stars with known properties through our pipeline for posterior
distributions, and generate Gaussian least squares fits to the output age
distributions. We generate these model stars at fixed Solar metallicity,
chosen at random with uniformly distributed ages between 1 and 50 Myr, and
masses from 0.1 to 1.3 M⊙. This mass limit roughly corresponds to where stars
at 1 Myr stop being clearly identifiable on the pre-main sequence as they
cross the subgiant branch. The photometry of these test stars is allowed to
vary according to the mean Gaia photometric uncertainty and mean reddening
uncertainty in the group, providing the expected photometric variation for an
individual star in that group. As for the Gaia sources, Gaussian distributions
are fit to the resulting output age and mass distributions for each test star,
and the input ages are binned according to the output fits for mass and age.
For very young stars where the peak is in the youngest age bin, a Gaussian
least squares fit will often view the probability distribution as the tail of
an exponential with a peak below zero, and therefore we exclude the 1 Myr bin
which often creates these anomalous fits, limiting the ages considered to
between 2 and 50 Myr. We compute a median input age and standard deviation for
the true ages of stars in each fit age and mass bin, which completes a map
that links any given set of age and mass results from the main pipeline to a
revised age result excluding any skewing effects related to metallicity and
multiplicity.

The vast majority of corrected age solutions are younger than those derived
from the peak of the Gaussian age distribution by less than $\sim$30-40%,
which is expected given that the addition of binaries adds a range of older
possible solutions to the photometry of the target. Some of the older (age
$>$30 Myr) and more massive stars near the limit of our method’s sensitivity
occasionally require larger correction to age of up to a factor of 2. A total
of 25727 of the original 30518 credible young candidates generate quality age
fits, or 84%. Most of the stars without good fits have photometric ages below
2 Myr, a category that includes both genuinely young stars and objects lowered
into that photometric age bracket by unresolved binarity or insufficiently
corrected reddening. We also lose some objects near the limits for
detectability where stars begin to merge into the main sequence. Upon
comparing the isochrones associated with these new revised ages to the
photometry of these stars, that the age results are consistent with our model.
Within Sco-Cen, many of the refined ages within subclusters host uncertainties
comparable to internal variation of the age between members, which would be
expected assuming a common formation time.

. The interference from binaries can however be reduced by considering groups
of stars rather than individuals. Therefore, when we search for age gradients
within a group, rather than considering the age fits of all stars
individually, we consider the median of the 10 nearest neighbors in space-
velocity coordinates. This is not a perfect solution for a group of stars, as
(see Section 3.5), however it is nonetheless a useful and efficient way to
identify and visualize age patterns in a population.

###  3.4 Clustering with HDBSCAN

The presence of structures in our young stellar population is clear through
the overdensities visible in Figure 7, however a consistent means of structure
identification is necessary for a proper analysis. To this end, we employ
HDBSCAN (Hierarchical Density-Based Spatial Clustering of Applications with
Noise; McInnes et al. 2017), a hierarchical density-based clustering algorithm
developed as an extension of the frequently-used DBSCAN clustering algorithm
(Ester et al., 1996). HDBSCAN has already been used quite successfully in
literature for a variety of purposes, including in Kounkel & Covey (2019),
where it is used to identify tenuous stellar structures of all ages in their
own Gaia DR2 sample. As such, there is strong precedent for its use in
identifying groups and associations within large stellar populations. Unlike
other clustering algorithms, which typically presuppose a density threshhold,
shape, or some other expected cluster property, HDBSCAN uses only a minimum
cluster size and a parameter with a smoothing effect on the density
distribution. As such, it is an algorithm with next to no built-in
assumptions, allowing for young structures to be identified without
discrimination based on complicating factors such as unusual shapes or wildly
varying scales. Furthermore, because the algorithm searches for structures at
different scales, it also provides a unique opportunity to explore different
levels of substructure within an association. Lastly, HDBSCAN enables
clustering in an arbitrary number of dimensions, allowing for both XYZ spatial
coordinates and l and b transverse velocities to be included in clustering
simultaneously. In some locations with high localized reddening, the Lallement
et al. (2019) maps are insufficient to properly correct for reddening,
resulting in anomalous overdensities without common velocities. Through the
inclusion of velocities in our clustering, we are able to exclude such
anomalies.

HDBSCAN’s search for overdensities makes use of a parameterized density
metric, which is based on the distance to the $k$th nearest neighbor (referred
to as the core distance, or $core_{k}($object$)$) (McInnes et al., 2017),
where $k$ is the first of two key HDBSCAN input parameters666given by the
HDBSCAN clustering parameter min_samples. Increasing $k$ effectively smooths
the density distribution and causes smaller clumps of stars to be ignored,
therefore serving as a proxy for how conservative the clustering will be.
HDBSCAN generates a new mutual reachability distance metric from these core
distances, which for objects $a$ and $b$ is defined as
$d_{mreach-k}(a,b)=max\\{core_{k}(a),core_{k}(b),d(a,b)\\}$, where $d(a,b)$ is
the distance from object $a$ to object $b$. Using this metric, HDBSCAN
generates a minimal spanning tree, which is essentially a web linking each
point in the data set through their nearest neighbor according to their mutual
reachability distances. By removing links in the tree of sources with the
lowest weights, corresponding to large mutual reachability distance to their
neighbors, links connecting clusters in the web are broken, leaving behind a
hierarchy of clusters that appear and fragment as the weight threshhold is
raised777An excellent visualization of this process is given on HDBSCAN’s
website, at https://hdbscan.readthedocs.io/en/latest/how_hdbscan_works.html.

The methods employed by HDBSCAN are comparable to those of DBSCAN, which
defines the cores of clusters by the presence of more than some number of
neighbors within a radius $\epsilon$, effectively taking only clusters that
exist at some arbitrary chosen scale (Ester et al., 1996). In HDBSCAN, rather
than looking for clusters at a single $\epsilon$, clustering is simultaneously
extracted at all scales, allowing the most persistent clusters to be
identified, rather than those that happen to emerge at any requested scale
(McInnes et al., 2017). Any clusters that come out of the tree with less than
$N$ members, where $N$ is the second key HDBSCAN input parameter888given by
the HDBCSAN clustering parameter min_cluster_size, are thrown out and merged
into their parent cluster: a larger-scale node on the clustering tree. The
quality of clusters in the tree is set by their persistence, which relates to
the range of scales over which the cluster is defined in the clustering tree.
HDBSCAN also includes the value $\epsilon$ as an optional parameter999given by
the HDBSCAN clustering parameter cluster_selection_epsilon, setting a minimum
threshold past which clusters can no longer be fragmented. This is useful in
cases where large groups with significant substructure are present, as it
allows for the entire region to be identified as a whole, without excessive
fragmentation. HDBSCAN has two options for cluster selection - “EOM” (excess
of mass), which selects the best or most persistent clusters in the clustering
tree, and “leaf”, which identifies only the nodes of the clustering tree,
effectively selecting clusters from a maximally fragmented clustering solution
(McInnes et al., 2017). EOM clustering is preferred in most cases, as the
persistence parameter that defines an EOM cluster is used by HDBSCAN as a
proxy for cluster quality, although leaf clustering can be useful when
attempting to detect subclusters within a well-defined higher-level cluster. A
more complete look at the methodology of HDBSCAN can be found in McInnes et
al. (2017), along with excellent visualizations of the tree clustering it
implements.

An additional consideration is the incorporation of velocities in clustering.
Since transverse velocity in km s-1 and galactic XYZ coordinates in pc do not
have matching scales, a corrective factor $c$ is required to

$(\text{X},\text{Y},\text{Z},c*v_{T,l},c*v_{T,b})$ (8)

The known young associations often have wildly different ratios between
spatial extent and velocity dispersions, however ratios of 4 to 6 pc/km s-1
are typical, approximately reflecting the ratios predicted for cloud scales of
25 to 100 pc in the Larson’s Law relation (Larson, 1981). While Larson’s Law
was designed to relate velocity dispersions to scales in clouds, associations
are expected to inherit properties from their parent cloud, and maintain those
properties for some time after formation, especially in loose associations
where dynamical interactions are infrequent (Larson, 1979; Wright & Mamajek,
2018; Ha et al., 2021). Size-velocity ratios of 4-6 pc/km s-1 predict
properties consistent with known groups such as Upper Sco and to a lesser
extent the more dispersed regions of UCL and LCC (Wright & Mamajek, 2018).
This range of factors is also consistent with relative sizes of visually-
identifiable Gaia clumps in spatial coordinates and in transverse velocity
space. Any variation within this range can be interpreted as a choice to tweak
the relative weight of velocity and spatial coordinates in the clustering,
where a larger corrective factor weighs velocity more heavily by reducing the
relative variability of the spatial component.

For all-sky clustering, a size-velocity corrective factor of $c$ = 6 pc/km s-1
was chosen, which weights kinematics slightly more heavily, reflecting the
larger velocity dispersions expected in the larger groups such as Sco-Cen that
we wish to identify at this level. A relatively large $\epsilon$ value of 25,
in units of pc (or $c$ km s-1 for velocity), supplements this choice,
preventing fragmentation of groups below a scale comparable to a moderately-
sized association like Upper Sco (Wright & Mamajek, 2018). This choice ensures
that groups like Sco-Cen can be treated like coherent singular groups, while
also allowing stellar structures that link subtly separated groups to be
detected, which would otherwise be incorporated into the background if those
groups were recorded separately. While groups being linked by this condition
does not guarantee that they share common formation origins, it does expose
spatial and kinematic similarities that can be further explored in
subclustering (see Section 5.4). The input parameter $N$, representing minimum
cluster size, was set to 10, as was $k$, which is set equal to $N$ by default.
This choice avoids the inclusion of negligibly small clusters or high-order
multiple systems, producing a set of largely visually convincing clusters.
However, some groups with less internally consistent velocities or spatial
distributions are also included, so we set a cluster persistence cut at a
value at 0.015. This number was chosen to ensure that no groups with velocity
distributions comparable to the field were included, most notably the
reddening anomaly towards Aquila marked in Figure 7, which is dense enough in
space to be marked as a cluster without any persistence-based quality
restrictions. Through this choice, our clustering results in this
implementation are made to be intentionally conservative, with the objective
of highlighting noteworthy structures in the Solar neighborhood.

We applied HDBSCAN using the above input parameters to members of our
population of probable young stars $\pi$/$\sigma_{\pi}>$25 A total of 28340
stars entered into the final 5-d data setHDBSCAN identified 27 spatially and
kinematically distinct young associations from this population, which are
shown in XY galactic spatial coordinates in Figure 10. These clustered sources
account for 41% of the total young population we identified, with the
remaining sources likely being composed of a mixture of young stars ejected
from larger associations, young stars in tenuous unrecognized associations,
and false positives enabled by reddening anomalies and unresolved
multiplicity. Sco-Cen was easily detected, representing by far the largest
group with just under 7400 stars, nearly an order of magnitude larger than the
next largest group, the near edge of Orion. Sco-Cen appears alongside a series
of nearby but spatially distinct structures, such as the dense clusters in
Chamaeleon. Perseus, Taurus, and the near edge of Orion are also identified in
at least some form, with varying levels of internal division, along with some
smaller known groups such as Perseus OB3 and the near edge of Vela. Many of
the remaining groups appear to be previously unknown.

Figure 10: The results of the all-sky clustering analysis from HDBSCAN in X/Y
galactic spatial coordinates (upper panel) and l/b sky coordinates (lower
panel). Each color corresponds to a separate top-level group identified by
HDBSCAN, which are labelled according to the IDs and names in Table 3. Groups
with substructure are shown using icons, which are given complete
subclustering analyses in Sections 4.1 (Taurus), 4.2 (Orion), 4.3 (Perseus),
4.4 (Sco-Cen and Chamaeleon), and 4.5 (all others). See Figure 7 for the X/Y
distribution of stars prior to clustering. Figure 11: The same as Figure 10,
but for transverse velocity.

To facilitate a more detailed study of Sco-Cen and other top-level groups we
performed an additional round of clustering on groups that showed visible
substructure, however the use of HDBSCAN there was modified slightly to enable
the recognition of hierarchical structures in each region. Rather than
identifying only the EOM clusters, as was done for top-level clustering, we
also extract the smaller-scale leaf clusters, which might exist within a
larger EOM cluster. The approach we employ fuses the two clustering methods,
recognizing EOM clusters containing leaf clusters as intermediate-level
groups, and the leaf clusters as subclusters within that larger EOM group. For
this clustering, $N$ and $k$ were kept the same, while the space-velocity
corrective factor $c$ was reduced to 5 to reflect the smaller scales, and
$\epsilon$ was removed entirely, allowing for subclusters to be identified as
separate regardless of the presence of nearby groups. The goal of this lower-
level clustering is to identify all potential subgroups in a structure,
regardless of scale or membership as part of a larger group. Ten of our 27
top-level subgroups have substructure, and within those ten groups this
implementation identifies a total of 60 EOM clusters, with eight of those 60
further subdividing into a total of 27 additional smaller leaf clusters. We
present the top level, EOM, and leaf clustering results for each candidate
young star in Table 1, alongside $P_{Age<50Myr}$ and the age solutions derived
in Section 3.3.

Table 1: The full catalog of all Gaia objects we identify with
$P(Age<50Myr)>0.1$, with basic Gaia properties, $P(Age<50Myr)$, and age as
derived in Section 3.3 for each. The TLC, EOM, and LEAF columns mark the top
level cluster, EOM cluster, and leaf cluster respectively that each star is a
part of, if any.

ID | TLCaaObjects with $\pi$/$\sigma_{\pi}<25$ were not subjected to clustering, and have a TLC and EOM markers of 0. Objects not assigned to a cluster are given the marker -1.  | EOMaaObjects with $\pi$/$\sigma_{\pi}<25$ were not subjected to clustering, and have a TLC and EOM markers of 0. Objects not assigned to a cluster are given the marker -1.  | LEAFbbleft blank if not related to a leaf cluster. | Gaia ID | RA (deg) | Dec (deg) |  $\pi$ (mas) | g | bp-rp | $P(Age<50Myr)$ | Ageccleft blank if no solution reached. (Myr)
---|---|---|---|---|---|---|---|---|---|---|---
|  |  |  |  |  |  |  |  |  |  | val | \+ err | \- err
0 | -1 | -1 |  | 3239529836638867840 | 73.8034 | 4.7175 | 4.170 | 17.80 | 3.50 | 0.34 | 28.3 | 8.4 | 5.7
1 | 22 | 11 |  | 6099439604319602560 | 217.3498 | -44.7853 | 5.960 | 15.60 | 3.00 | 0.94 | 12.0 | 2.7 | 3.0
2 | -1 | -1 |  | 4154795132763964032 | 278.9114 | -10.8489 | 3.020 | 14.40 | 2.40 | 1.00 |  |  |
3 | 21 | 4 |  | 5774202930948040064 | 266.3687 | -82.1980 | 6.260 | 16.10 | 3.10 | 0.90 | 18.2 | 5.6 | 4.8
4 | -1 | -1 |  | 3442403338420195840 | 82.6082 | 28.2391 | 3.460 | 11.90 | 1.40 | 0.23 |  |  |
5 | 17 | -1 |  | 3442418800302357248 | 83.0898 | 28.3245 | 5.220 | 16.10 | 3.00 | 0.89 | 15.6 | 3.3 | 5.3
6 | -1 | -1 |  | 4154777227001098112 | 279.5766 | -10.6085 | 3.100 | 14.30 | 1.90 | 0.15 |  |  |
7 | -1 | -1 |  | 1937498784187322624 | 351.8766 | 45.0985 | 3.830 | 16.70 | 2.90 | 0.24 | 22.8 | 9.6 | 4.3
8 | -1 | -1 |  | 6195524512419762816 | 204.7577 | -21.6912 | 11.930 | 13.30 | 2.60 | 0.56 |  |  |
9 | -1 | -1 |  | 1716855559591213824 | 192.1970 | 78.9078 | 7.670 | 15.10 | 2.90 | 0.46 | 19.9 | 4.8 | 4.8

††footnotetext: only a small subset of the 30518 total objects is shown here.
The full table is available in the online version of this paper.

###  3.5 Group Ages

The ages we derive for entire groups of stars cannot be calculated in the same
way as was done for individual stars, as the population we consider during
clustering is, by design, skewed towards photometrically younger objects. In
groups of fixed age, photometrically younger stars include objects with
unresolved companions, or more generally stars on the high side of the
expected photometric variation. This means that group ages calculated using
our likely young stars alone will be skewed young, especially in groups older
than $\sim$20 Myr . We can reduce this effect by reintroducing stars that our
pipeline did not identify as young, and fitting an isochrone to the expanded
population. This can be done by searching for objects that have spatial
coordinates and kinematics consistent with the known members of that group, at
least as close to the tenth-nearest HDBSCAN-identified member in space-
velocity coordinates as the most peripheral HDBSCAN-identified member,
mimicking the original HDBSCAN methods for finding members of clusters with
$k$=10. This approach effectively uses the groups in the original clustering
analysis as signposts from which more complete young stellar populations can
be gathered and used for age estimation, removing influences from our
pipeline’s selection biases. The stars present in these extended populations
are provided in Table 2.

For each candidate young cluster, we perform a single age fit based on the
photometry of , including both the original candidate young stars and the
extended population of comoving/cospatial objects mentioned above. Our methods
employ a least squares fit, where we compare the photometry of these
populations to a grid of solar metallicity PARSEC isochrones with ages between
0.25 and 80 Myr (Chen et al., 2015). Only stars within the range of G
magnitudes occupied by the isochrone grid are included to enable reliable
interpolation, a restriction that also helps to suppress some of the
background contamination added upon reintroducing candidate members that are
not photometrically young. We also to restrict the sample to the pre-main
sequence. Most of the groups we recover have pre-main sequences that fall
comfortably within the range of photometry given by these 0.25-80 Myr
isochrones, with the exception of the Pleiades, which we exclude from age
estimation due to its considerably older known age (e.g., Lodieu et al.,
2019).

To reduce contamination from external sources, we restrict the sources
included in our age fits according to the core distance (as defined in Section
3.4). For the core distance restrictions, we define a weight parameter in
which the group member with the smallest core distance in the association has
a weight of 1, the member with the largest core distance has a weight of 0,
and sources in between have weights spaced linearly according to the core
distance between these extremes. For top-level clustering, a default condition
of weight $>$ 0.3 was required for the inclusion of members in the age fit for
the group, which in most cases ensures that a young pre-main sequence
dominates background contamination. In more contaminated groups this was
raised as high as 0.5. For subclusters, we often found that no weight cut was
required due to tighter spatial and kinematic distributions that are less
vulnerable to contamination, however we re-imposed the restrictions where
necessary. Very few groups needed changes to these conditions, with some
examples being the large and badly contaminated subgroups of CHA-1 and ORI-2,
as well as Taurus-Orion IV and Cepheus Flare at the top level.

further restricted the population used for age fitting based on Re-normalized
Unit Weight Error (RUWE), (e.g., Bryson et al., 2020). In most cases, we
removed sources with RUWE$>$1.1, roughly corresponding to where the additional
contribution to the RUWE curve from binaries becomes negligible in Bryson et
al. (2020). This harsh RUWE cut works well to remove binary sequences, however
in cases with consistently high RUWE, such groups with embedded sources (e.g.
Perseus, see Sec. 4.3), this condition occasionally removes nearly all stars.
Therefore, in cases where more than two-thirds of stars are removed by
RUWE$<$1.1, we relax the condition to RUWE$<$1.2. This value corresponds to
approximately where Gaia populations are expected to become fully binary-
dominated, and it is the recommended criterion for culling binaries in Bryson
et al. (2020). This cut combined with the weight cut significantly reduced the
number of new stars used in the fits, from a typically order-of-magnitude
increase in populations upon the initial reintroduction of objects to a number
more consistent with the expected number, which should be equal to the number
of young candidates times the inverse of the expected recovery rate (see
Section 3.2).

Once these culling conditions were applied, visually convincing fits were
produced for all groups in our sample, which neglect strong background
contamination and visible binary sequences. Figure 12 shows a two examples of
the results of this fitting. We generate fits for all HDBSCAN-identified
groups without substructure, as well as a few cases that do show substructure
but with a limited enough age spread for a global age to make sense. As we
note throughout Section 4, the vast majority of our age fits for known groups
are consistent with literature values, however the uncertainties should be
treated with caution, as we report standard errors that cannot capture the
dominant systematic errors related to the choice of model and clipping, or any
inherent age spread. Due to the 80 Myr isochrone limit for inclusion in
fitting, some members in associations with ages $>$50 Myr where these
isochrones become close together may be removed, so some additional caution
should be applied to the age fits near the upper end of our age range.

Table 2: The full catalog of Gaia stars with space-velocity positions
consistent with our photometrically young populations. The TLC, EOM, and LEAF
columns mark the Top Level cluster, EOM cluster, and LEAF cluster that each
star is a part of, if any. We also include basic Gaia properties and the
weight parameter described in Section 3.5, which serves as a proxy for
closeness of the group center.

IDaaNote that in complex environments, extreme outliers may be shared between two clusters. These cases receive one entry for each possible solution. | TLC | EOMbbObjects not assigned to an EOM cluster are given the marker -1.  | LEAFccLeft blank if not related to a leaf cluster. | Gaia ID | RA | Dec | parallax | g | bp-rp | weight
---|---|---|---|---|---|---|---|---|---|---
|  |  |  |  | (deg) | (mas) | (mag) |
0 | 1 | -1 |  | 2297218611804127104 | 309.8679 | 79.8920 | 3.176 | 18.93 | 2.87 | 0.02
1 | 1 | -1 |  | 2226243658763815552 | 336.1287 | 69.4314 | 3.035 | 17.39 | 3.61 | 0.23
2 | 1 | -1 |  | 2226228579136886144 | 336.4913 | 69.3191 | 3.296 | 13.51 | 1.23 | 0.08
3 | 1 | -1 |  | 2286800434190203776 | 334.4768 | 80.5854 | 3.191 | 17.79 | 2.69 | 0.09
4 | 1 | -1 |  | 2297464666187525120 | 324.2971 | 80.7429 | 3.391 | 11.98 | 0.83 | 0.01
5 | 1 | -1 |  | 2212189293119019264 | 343.6841 | 66.3300 | 3.134 | 17.04 | 2.58 | 0.02
6 | 1 | -1 |  | 2226302723153297664 | 336.9986 | 69.7687 | 3.066 | 18.36 | 3.76 | 0.27
7 | 1 | -1 |  | 2286846781182031744 | 339.9039 | 80.6398 | 3.300 | 15.34 | 1.98 | 0.08
8 | 1 | -1 |  | 2229484297488789248 | 338.5122 | 70.9616 | 3.439 | 17.85 | 2.65 | 0.05
9 | 1 | -1 |  | 2229446020739020032 | 336.3050 | 70.5487 | 3.399 | 17.84 | 2.56 | 0.05

††footnotetext: only a small subset of the 94424 total objects is shown here.
The full table is available in the online version of this paper.

Figure 12: Sample isochrone fits for two of our young associations. In the
left panel, we display our handling of severe contamination in Taurus-Orion
IV, where a weight cut raised to 0.5 resulted in a very clean pre-main
sequence emerging from the noise. Objects that failed the weight cuts are
marked with red circles, and the stars that pass are given by black diamonds.
In the right panel, we present members of the relatively uncontaminated
Carina-Musca association after applying weight cuts. Stars marked by large
cyan circles fail the subsequent RUWE cut, while stars marked by black
diamonds pass it. Stars that fail the RUWE cut are likely binaries, hence the
existence of a major concentration of them along an apparent binary sequence.
Only objects marked by orange circles are identified as young Carina-Musca
members using our methods, highlighting the potential bias towards
photometrically younger stars in that population that our age fitting method
attempts to rectify.

##  4 Results

By applying HDBSCAN clustering to our sample of young stars, we have
identified extensive young structures in the solar neighborhood. These include
27 top-level regions (see Fig. 10), with properties provided in Table 3. Many
of these regions have significant internal structure, often with well-
separated subclusters, or subclusters widely separated in age. Performing a
subclustering analysis on each of these regions allows for insight into their
star formation histories. We therefore manually inspect each region for
visible substructure, and apply HDBSCAN subclustering on the ten regions where
substructure is visually evident. While we do not assert an absence of
subclustering in the remaining regions, most are too small or sparse for a
subclustering analysis to be useful, and will therefore likely require a more
targeted census before substructure can be ascertained. Details on this
HDBSCAN implementation are described in Section 3.4. Further detail on
individual complex regions will be left to future publications, as will
precise characterization of some of the previously unknown or little-known
groups that we identify for the first time in this publication. These lesser-
known groups are discussed in Section 4.5, along with a discussion of their
relations to past literature. Most are assigned a unique name in Table 3 to
reflect their distinctiveness as significant nearby young associations.

Table 3: Properties of the nearby clusters identified using HDBSCAN. NC in the
age column indicates that we find the group to be non-coeval, see the age
column in the relevant subclustering table for the ages of the subgroups
(e.g., Tables 4, 5, 6, and 7).

TLC | Name | N | RA | Dec | l | b |  $D_{sky}$aaOn-sky spatial extent in galactic l/b, in the form of the RMS in major axis $\times$ minor axis, when fit with a bivariate Gaussian. | d | $\mu_{RA}$ | $\mu_{Dec}$ | $V_{T,l}$ | $V_{T,b}$ |  $\sigma_{V_{T}}$bbThe radial extent of the velocity distribution, in the form of the RMS in semi-major axis $\times$ semi-minor axis, when fit with a bivariate Gaussian. | Age
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
|  |  | (deg) | (deg) | (deg) | (pc) | (mas/yr) | (km/s) | (km/s) | (Myr)
1 | Cepheus Flare | 17 | 339.5 | 73.9 | 113.9 | 13.4 | 5.3$\times$2.2 | 329.3$\pm$3.2cchas a member within 1 pc of the search horizon, and therefore members are likely present beyond that limit. | 6.8 | 0.9 | 9.9 | -4.0 | 2.0$\times$1.4 | 11.2$\pm$3.3
2 | PleiadeseeOlder cluster - Binary sequence dominant in detection. | 23 | 57.4 | 24.9 | 166.4 | -22.4 | 7.6$\times$4.4 | 139.9$\pm$6.7 | 18.9 | -44.7 | 29.2 | -13.4 | 1.3$\times$1.2 |  $>$80
3 | Taurus-Orion I$\dagger$$\dagger$footnotemark:  | 24 | 65.9 | 13.5 | 181.6 | -24.4 | 5.3$\times$2.5 | 299.2$\pm$10.7 | 4.5 | -5.2 | 9.7 | 0.4 | 2.0$\times$0.6 | 10.6$\pm$2.2
4 | Ophiuchus Southeast$\dagger$$\dagger$footnotemark:  | 31 | 257.5 | -18.3 | 4.5 | 12.7 | 4.7$\times$1.7 | 217.5$\pm$14.2 | -5.4 | -11.6 | -13.0 | -2.2 | 1.6$\times$0.7 | 17.5$\pm$1.9
5 | Fornax-Horologium$\dagger$$\dagger$footnotemark:  | 40 | 54.2 | -34.5 | 235.1 | -53.7 | 20.1$\times$9.6 | 107.8$\pm$10.7 | 34.7 | -3.9 | 1.9 | 17.3 | 3.9$\times$1.7 | 38.5$\pm$3.6
6 | CMa North$\dagger$$\dagger$footnotemark:  | 28 | 102.3 | -15.3 | 226.4 | -7.4 | 5.5$\times$4.8 | 184.9$\pm$19.7 | -0.2 | -5.4 | 4.1 | -2.3 | 1.7$\times$0.8 | 30.2$\pm$2.7
7 | Aquila East$\dagger$$\dagger$footnotemark:  | 30 | 297.8 | -9.0 | 31.6 | -17.3 | 12.2$\times$3.8 | 136.6$\pm$6.5 | 8.5 | -25.9 | -12.7 | -12.2 | 1.9$\times$0.8 | 20.2$\pm$1.5
8 | Cepheus Far North$\dagger$$\dagger$footnotemark:  | 219 | 329.6 | 74.7 | 111.7 | 17.8 | 16.9$\times$5.7 | 179.1$\pm$26.3 | 14.8 | 6.1 | 14.1 | -2.9 | 2.4$\times$0.9 | 23.8$\pm$1.4ddHas substructure, but approximately coeval.
9 | Vela-CG7 | 111 | 124.9 | -43.7 | 260.8 | -4.2 | 7.3$\times$4.6 | 323.7$\pm$9.0cchas a member within 1 pc of the search horizon, and therefore members are likely present beyond that limit. | -7.3 | 10.1 | -19.2 | -0.5 | 1.2$\times$0.9 | 14.4$\pm$1.2
10 | ASCC 123 | 34 | 334.0 | 53.8 | 101.3 | -2.2 | 8.0$\times$3.5 | 251.1$\pm$24.4 | 12.0 | -2.0 | 10.3 | -9.8 | 2.2$\times$1.1 | 49.5$\pm$2.2
11 | Cepheus-Cygnus$\dagger$$\dagger$footnotemark:  | 36 | 311.7 | 56.6 | 93.8 | 8.4 | 5.9$\times$2.6 | 314.2$\pm$11.7cchas a member within 1 pc of the search horizon, and therefore members are likely present beyond that limit. | 5.1 | 3.6 | 8.9 | -2.7 | 1.0$\times$0.7 | 32.6$\pm$1.2
12 | Lyra$\dagger$$\dagger$footnotemark:  | 29 | 282.7 | 32.8 | 62.7 | 14.4 | 3.5$\times$1.9 | 300.1$\pm$11.9 | 1.7 | -3.9 | -4.1 | -4.4 | 1.0$\times$0.6 | 31.0$\pm$1.8
13 | Cerberus$\dagger$$\dagger$footnotemark:  | 66 | 271.4 | 22.1 | 48.3 | 19.5 | 9.3$\times$2.9 | 308.6$\pm$13.3cchas a member within 1 pc of the search horizon, and therefore members are likely present beyond that limit. | -1.0 | -7.3 | -10.4 | -2.7 | 2.2$\times$0.9 | 30.0$\pm$0.8
14 | Carina-Musca$\dagger$$\dagger$footnotemark:  | 168 | 158.8 | -69.4 | 291.0 | -10.8 | 15.0$\times$3.3 | 241.8$\pm$33.4 | -25.8 | 3.1 | -27.7 | -12.5 | 1.1$\times$0.9 | 32.0$\pm$1.0ddHas substructure, but approximately coeval.
15 | Perseus | 264 | 58.4 | 32.9 | 161.4 | -15.9 | 6.4$\times$4.1 | 294.8$\pm$22.4cchas a member within 1 pc of the search horizon, and therefore members are likely present beyond that limit. | 5.3 | -8.1 | 13.0 | -3.6 | 2.5$\times$1.2 | NC
16 | Perseus OB3 | 145 | 56.1 | 44.1 | 153.2 | -7.7 | 21.2$\times$6.7 | 167.5$\pm$18.6 | 22.5 | -26.6 | 26.7 | -5.8 | 2.4$\times$0.9 | NC
17 | Taurus-Orion II$\dagger$$\dagger$footnotemark:  | 64 | 89.5 | 18.3 | 190.6 | -2.8 | 13.1$\times$4.1 | 207.4$\pm$17.2 | 7.6 | -18.3 | 19.1 | -2.6 | 2.5$\times$0.8 | 41.5$\pm$1.7
18 | Greater Taurus | 674 | 73.3 | 19.1 | 181.6 | -14.9 | 23.9$\times$13.5 | 142.8$\pm$24.9 | 7.6 | -22.1 | 14.5 | -4.7 | 5.0$\times$3.4 | NC
19 | IC 2391 | 101 | 139.1 | -58.0 | 276.7 | -8.3 | 23.2$\times$4.2 | 157.7$\pm$15.8 | -26.7 | 18.9 | -24.8 | -4.3 | 1.8$\times$0.7 | 45.0$\pm$1.7
20 | NGC 2451A | 51 | 116.7 | -40.9 | 255.1 | -8.0 | 10.0$\times$2.4 | 189.0$\pm$14.0 | -22.2 | 16.4 | -22.3 | -10.1 | 1.4$\times$0.9 | 48.5$\pm$1.3
21 | Chamaeleon | 259 | 188.0 | -74.3 | 301.8 | -13.0 | 14.9$\times$11.6 | 190.1$\pm$14.5 | -18.2 | -5.0 | -17.7 | -7.7 | 3.3$\times$1.2 | NC
22 | Sco-Cen | 7394 | 223.0 | -40.5 | 329.1 | 9.4 | 50.4$\times$20.0 | 137.8$\pm$24.3 | -18.5 | -19.3 | -18.0 | -6.0 | 2.6$\times$2.0 | NC
23 | Taurus-Orion III$\dagger$$\dagger$footnotemark:  | 27 | 76.4 | 11.8 | 189.4 | -17.2 | 5.0$\times$2.5 | 309.6$\pm$10.1 | 1.9 | -3.5 | 5.9 | -0.5 | 0.7$\times$0.6 | 15.2$\pm$2.8
24 | Vela-CG4 | 551 | 108.7 | -41.3 | 252.8 | -13.9 | 21.5$\times$8.8 | 265.7$\pm$45.0cchas a member within 1 pc of the search horizon, and therefore members are likely present beyond that limit. | -11.0 | 7.5 | -14.1 | -8.7 | 2.9$\times$0.9 | 33.7$\pm$0.9ddHas substructure, but approximately coeval.
25 | Taurus-Orion IV$\dagger$$\dagger$footnotemark:  | 108 | 73.1 | 19.4 | 181.0 | -15.0 | 18.1$\times$8.7 | 218.8$\pm$27.4 | 1.4 | -7.6 | 7.1 | -3.9 | 2.8$\times$1.6 | 21.3$\pm$2.1
26 | Monoceros Southwest$\dagger$$\dagger$footnotemark:  | 292 | 100.3 | -4.6 | 215.9 | -4.4 | 8.5$\times$4.1 | 279.0$\pm$35.2cchas a member within 1 pc of the search horizon, and therefore members are likely present beyond that limit. | -6.6 | -2.7 | -0.7 | -9.3 | 1.7$\times$0.9 | 25.5$\pm$0.6ddHas substructure, but approximately coeval.
27 | Greater Orion | 757 | 79.7 | -5.5 | 207.1 | -23.0 | 12.5$\times$9.1 | 299.1$\pm$34.0cchas a member within 1 pc of the search horizon, and therefore members are likely present beyond that limit. | 1.0 | -1.3 | 2.1 | 0.4 | 1.8$\times$1.3 | NC

††footnotetext: Name is new - limited to no coverage in the literature or left
unnamed in the literature that does exist.

###  4.1 Greater Taurus

Table 4: Clustering within the Greater Taurus. GT is the Greater Taurus ID, a number from EOM clustering used for identification. Leaf clusters within an EOM cluster are given a subgroup ID, which is given in the LEAF column. NC in the age column marks subclustered groups with non-coeval populations. GT | LEAF | Name | N | RA | Dec | l | b |  $D_{sky}$aaOn-sky spatial extent in galactic l/b, in the form of the RMS in major axis $\times$ minor axis, when fit with a bivariate Gaussian. | d | $\mu_{RA}$ | $\mu_{Dec}$ | $V_{T,l}$ | $V_{T,b}$ |  $\sigma_{V_{T}}$bbThe radial extent of the velocity distribution, in the form of the RMS in semi-major axis $\times$ semi-minor axis, when fit with a bivariate Gaussian. | Age
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
|  |  |  | (deg) | (deg) | (deg) | (pc) | (mas/yr) | (km/s) | (km/s) | (Myr)
1 |  |  $\mu$ Tau | 17 | 59.0 | 9.6 | 180.0 | -32.0 | 8.5$\times$2.8 | 153.5$\pm$6.1 | 23.7 | -23.3 | 24.1 | 1.9 | 1.3$\times$1.0 | 52.2$\pm$3.4
2 |  |  | 25 | 80.3 | -6.2 | 208.1 | -22.9 | 6.0$\times$3.5 | 144.2$\pm$6.2 | -1.1 | -3.0 | 1.5 | -1.6 | 0.8$\times$0.6 | 31.0$\pm$2.5
3 |  |  | 46 | 77.6 | -2.0 | 202.6 | -23.2 | 11.1$\times$5.2 | 165.5$\pm$8.2 | 10.0 | -14.2 | 13.6 | 1.5 | 2.4$\times$1.3 | NC
3 | A |  | 13 | 78.1 | 2.0 | 199.1 | -20.9 | 5.2$\times$1.6 | 166.9$\pm$3.6 | 8.6 | -14.1 | 13.1 | 0.3 | 0.7$\times$0.6 | 19.3$\pm$3.9
3 | B |  | 12 | 80.8 | -8.4 | 210.5 | -23.4 | 2.0$\times$1.2 | 157.2$\pm$2.3 | 8.3 | -12.4 | 11.0 | 1.5 | 0.4$\times$0.4 | 13.8$\pm$2.2
4 |  |  | 51 | 68.0 | 19.7 | 177.6 | -18.8 | 11.2$\times$6.6 | 121.6$\pm$7.8 | 0.7 | -15.1 | 7.0 | -5.3 | 1.6$\times$1.0 | NC
4 | A |  | 16 | 70.6 | 20.4 | 178.7 | -16.6 | 7.0$\times$2.9 | 119.8$\pm$4.4 | -1.1 | -16.1 | 6.7 | -6.2 | 0.7$\times$0.5 | 37.8$\pm$3.2
4 | B |  | 11 | 61.1 | 20.1 | 172.7 | -23.6 | 3.3$\times$2.8 | 120.5$\pm$3.3 | 4.8 | -14.0 | 7.7 | -3.5 | 0.5$\times$0.3 | 29.8$\pm$8.3
5 |  | Theia 93 | 55 | 80.7 | 24.2 | 181.3 | -6.8 | 4.3$\times$3.1 | 177.3$\pm$7.4 | 1.8 | -18.1 | 13.5 | -7.2 | 0.9$\times$0.6 | 8.5$\pm$1.3
6 |  | 118 Tau | 33 | 84.2 | 23.0 | 184.0 | -4.8 | 5.8$\times$4.5 | 108.8$\pm$3.6 | 6.1 | -37.6 | 18.1 | -7.7 | 0.8$\times$0.6 | 16.8$\pm$2.0
7 |  | 32 Ori | 11 | 81.5 | 8.3 | 195.3 | -14.7 | 3.6$\times$3.3 | 102.0$\pm$5.1 | 8.8 | -33.2 | 15.9 | -4.6 | 0.7$\times$0.6 | 27.2$\pm$3.8
8 |  | L1527ccThis named structure is not included in any HDBSCAN-defined leaf subclusters. Names associated with subgroups are also contained in the higher-level EOM group. | 87 | 66.6 | 26.4 | 171.5 | -15.5 | 5.5$\times$1.8 | 130.7$\pm$4.5 | 7.9 | -23.8 | 14.2 | -6.2 | 1.7$\times$0.7 | 4.5$\pm$0.9ddHas substructure, but approximately coeval.
8 | A | L1529 | 30 | 68.3 | 24.5 | 174.0 | -15.6 | 3.1$\times$1.5 | 129.4$\pm$2.7 | 7.3 | -21.5 | 12.9 | -5.1 | 0.7$\times$0.6 | 3.4$\pm$0.5
8 | B | B209/L1495 | 38 | 65.0 | 28.0 | 169.3 | -15.6 | 2.6$\times$0.7 | 131.2$\pm$3.3 | 8.6 | -25.6 | 15.2 | -7.0 | 0.7$\times$0.4 | 3.8$\pm$0.8
9 |  | L1544ccThis named structure is not included in any HDBSCAN-defined leaf subclusters. Names associated with subgroups are also contained in the higher-level EOM group. | 77 | 77.1 | 29.0 | 175.4 | -6.7 | 7.1$\times$3.7 | 159.2$\pm$5.6 | 4.4 | -25.3 | 17.4 | -8.6 | 1.0$\times$0.6 | NC
9 | A | L1517 | 40 | 75.5 | 30.1 | 173.7 | -7.1 | 4.1$\times$1.8 | 158.0$\pm$3.3 | 4.3 | -25.1 | 16.9 | -8.8 | 0.7$\times$0.4 | 5.6$\pm$0.8
9 | B |  | 10 | 81.4 | 25.8 | 180.2 | -5.4 | 1.6$\times$0.7 | 162.1$\pm$3.6 | 3.6 | -25.5 | 17.8 | -8.5 | 0.4$\times$0.2 | 11.5$\pm$1.3
10 |  | L1551 | 30 | 69.0 | 17.5 | 180.2 | -19.7 | 3.1$\times$1.8 | 145.6$\pm$3.8 | 12.3 | -19.1 | 15.6 | -1.5 | 0.7$\times$0.4 | 7.0$\pm$2.0
11 |  | B213/L1536 | 34 | 67.4 | 24.9 | 173.2 | -16.0 | 5.3$\times$1.1 | 159.2$\pm$4.2 | 10.7 | -17.4 | 15.2 | -2.5 | 0.8$\times$0.7 | 3.3$\pm$0.9
Figure 13: Members of the EOM groups that HDBSCAN identifies in Greater Taurus
in RA/Dec, galactic X/Y, and transverse velocity coordinates. icons indicate
known Taurus members and icons indicate objects not previously recognized as
part of Taurus. Small grey dots mark unclustered Taurus members. Each group is
labelled according to its Greater Taurus ID from Table 4 in the top panel. Our
identification of Greater Taurus includes significant populations beyond the
known membership. The main stars of the Taurus and Orion constellations are
included for reference in the top panel, with the dotted circle attached to
Taurus representing the Pleiades. Figure 14: Subclustering in the four Greater
Taurus EOM clusters where it is present, displayed in RA/Dec sky coordinates.
In all cases, unclustered members of the parent group are marked as small grey
dots, and the subgroups are given as squares and diamonds. The ID of the
parent EOM cluster is indicated on each panel.

The general direction towards the Taurus constellation is known to contain
multiple nearby stellar populations. The largest of these is associated with
sites of active star formation in the Taurus Molecular Cloud (TMC), which form
an extensive network of young stars often referred to as the Taurus-Auriga
Association or simply “Taurus”, as we refer to it here (e.g., see Kenyon et
al., 2008). This complex represents one of the nearest and therefore most
accessible sites of active star formation to the Sun ($\sim$140 pc, Galli et
al., 2018), and as such it has been subject to many recent studies to expand
the known stellar populations, learn more about the structure and history of
the region, and study star formation more broadly . Through some of that work,
kinematic substructure has recently been identified in the region, dividing
the stellar populations into four different subgroups, each tied to a
different set of clouds within the TMC (Luhman, 2018). The lists of known and
proposed Taurus members continue to expand, although the existence of older
members beyond the well-established young ($<$10 Myr) populations remains
controversial (Kraus et al., 2017; Luhman, 2018). Aside from the Taurus
populations associated with the TMC, there are also many smaller and less-
studied groups in the area, including the 118 Tauri, 32 Orionis, and $\mu$
Tauri Associations (Mamajek, 2016, 2007; Gagné et al., 2020). Since all of
these groups have only recently been fully characterized, any relations these
associations might have to one another and the other major stellar populations
in Taurus remain unsettled.

Using HDBSCAN, we identify an extensive contiguous network of young stars in
the direction of Taurus, forming a top-level group that we refer to as
“Greater Taurus”. Significant substructure is found within the region,
breaking it into a total of 11 EOM groups which are displayed in Figure 13.
The velocity dispersions of Taurus subgroups are typically of order 1 km s-1
or smaller, with a larger dispersion for the group as a whole, at about 4 km
s-1. However, most groups closely associated with the TMC are in a somewhat
tighter distribution, with even the most extreme exemplars differing in
transverse velocity by less than 10 km -1. Outlying groups typically differ
from the core TMC groups in transverse velocities by 5-15 km s-1, however
their distributions remain connected to each other and to the core Taurus
groups by lower density distributions of stars in the intervening parameter
space, hence HDBSCAN merges them into a single large-scale Greater Taurus
group (see Fig. 13). Since the area Taurus covers is quite large, projection
effects account for a significant fraction of the velocity changes. Assuming a
constant velocity equal to the Taurus UVW vector from Luhman et al. (2009)
would generate a 3 km s-1 velocity dispersion from projection alone. There
are, however, subgroups with notably different velocities compared to the
projected Taurus UVW vector, so internal velocity spreads do exist. Known
members in Taurus are marked using icons in Figure 13 as a way to compare our
sample to the currently known extent of the region (Kraus et al., 2017;
Luhman, 2018; Esplin & Luhman, 2019; Galli et al., 2019; Krolikowski,
submitted). The subgroups we identify as part of Greater Taurus extend far
beyond these known Taurus members. Higher-order subclusters are also found at
the leaf level, with 4 of 11 EOM groups dividing further into 2 leaf groups
each. Close-up views of the groups containing leaf subclusters are provided in
Figure 14.

Four of the groups that HDBSCAN identifies appear to be associated with the
ongoing and recent star formation in the TMC (e.g., Rebull et al., 2011;
Luhman, 2018). These groups contain the majority of currently known Taurus
Association members (see Figure 14), and they all correspond to known
populations of Young Stellar Objects (YSOs) in and around active dense clouds
in the TMC. GT-8 contains YSOs associated with L1495 B209, L1527, and L1529;
GT-9 contains L1517 and L1544; GT-10 contains L1551; and GT-11 contains B213
and L1536 (Luhman, 2018; Onishi et al., 2002). These four groups are similar
to those proposed in Luhman (2018), with the main differences being our
separation between GT-10 and GT-11, which Luhman (2018) groups together into
their “blue” cluster, and our exclusion of their “green” cluster as a distinct
subgroup due to its very small membership. All four of our core TMC groups
consist almost entirely of known Taurus member stars, with the notable
exception of a southeastern extension to GT-9, which forms its own leaf
cluster, GT-9B (recognized as Theia 68 in Kounkel & Covey 2019). The other
GT-9 subgroup, GT-9A, consists of the known YSOs in the region around L1517.
GT-8 also fragments into leaf clusters, with GT-8A consisting mainly of YSOs
associated with L1529, and GT-8B consisting of B209 and L1495. All of these
core TMC groups (aside from GT-9B) have derived ages between 3 and 7 Myr,
consistent with our knowledge that active star formation is happening in these
regions.

Distributed populations around Taurus have long been proposed, particularly in
the direction of GT-1 (e.g., Neuhaeuser et al., 1995; Wichmann et al., 1996;
Magazzu et al., 1997), however these early papers have reached mixed results
on potential association with the TMC (e.g., Briceno et al., 1997). We also
identify distributed populations in our work, forming an extensive network of
peripheral groups that we associate with Greater Taurus, many of which are
significantly older than the core TMC groups. Some of these groups, including
GT-5, GT-6, and both subgroups of GT-4 contain proposed Taurus members, but as
indicated in Figure 13, their coverage in surveys of Taurus to date is far
from complete (Krolikowski, submitted). The youngest of these groups is GT-5
(which has considerable shared membership with Theia 93 from Kounkel & Covey,
2019), a subgroup with TMC-consistent kinematics that appears to represent a
slight continuation of Taurus to farther distances. Others have ties to the
TMC that are considerably less clear. GT-6 has proposed Taurus members in it,
however it was previously identified independently as the 118 Tauri Moving
Group by Mamajek (2016), which has membership falling entirely within GT-6 and
kinematically consistent objects. This group represents one of three we
identify that have previously been identified independently from the core
populations in Taurus, with the other two being GT-1 and GT-7. In GT-7, over
half of our membership overlaps with that of the 32 Orionis Association
(Mamajek, 2007; Bell et al., 2017), a somewhat older group where our derived
age of 27.2 $\pm$ 3.8 Myr is in agreement with the 22$\pm$4 Myr isochronal age
solution from Bell et al. (2015). The remainder of our membership matches
spatially and kinematically with known 32 Orionis members, and while we only
identify 11 members of this group, its established membership from Bell et al.
(2017) contains only 47 objects, making the sample we gather smaller, but
consistent with recovery rates for groups in this age range (see Section 3.2).

Our member stars for GT-1 match with those of the much older $\mu$ Tauri
Association in both spatial coordinates and kinematics, where our derived age
of $\sim$52.2$\pm$3.4 Myr is marginally within uncertainties of the 62$\pm$7
Myr solution derived by Gagné et al. (2020). The Gagné et al. (2020) candidate
member list contains all of our GT-1 members. With only 17 GT-1 members
compared to $\sim$500 $\mu$ Tau candidates from Gagné et al. (2020), our
detection of this group is limited, although this low recovery rate is likely
primarily due to the very old age of the region (see Section 3.2). A few GT-3
members in its far western extension are also included in the Gagné et al.
(2020) catalog, an overlap indicative of a possible link between GT-1 ($\mu$
Tauri) and GT-3. The cores of these two groups are separated by $\sim$60 pc in
space, $\sim$10 km s-1 in transverse velocity, and $>$30 Myr in age, and
therefore their affiliation with one another is far from certain. However, the
separation of these subgroups in velocity is consistent with projection
effects, and a low-density, near-contiguous linear structure in the plane of
the sky does appear to connect them, lending plausibility to the connection
between these subgroups.

Aside from this slight overlap with $\mu$ Tau, the remainder of the stars in
GT-3 and GT-2 do not match with any current Taurus catalogs related to the
TMC, as expected given that these groups stretch up to 25 degrees south of the
core TMC groups to the crowded fields in front of the Orion Association. GT-2
(along with GT-4) is grouped under Kounkel & Covey (2019) group Theia 44, but
these groups are otherwise absent in literature. GT-3 divides into two
subgroups: GT-3A to the north, and GT-3B to the south, with ages of $\sim$19
and 14 Myr, respectively. GT-3B is projected directly behind the somewhat
older ($\sim$21 Myr) GT-2 group in the plane of the sky, however despite this
close physical proximity between GT-3B and GT-2, they are separated by more
than 10 km s-1 in transverse velocity space, making close association in
formation origin unlikely. GT-3, together with the GT-4 and GT-1, represent
visible outliers from the relatively compact transverse velocity distribution
occupied by the core TMC groups, which is less than 10 km s-1 across. While
GT-2 also qualifies as an outlier to the TMC in transverse velocities, it is
actually not inconsistent with a common origin there, as its velocities are
consistent with the projection of the Taurus UVW velocity vector at its
location (Luhman et al., 2009). GT-4 is unique in that it has an outlying
kinematic signature despite it physically overlapping with GT-8 and GT-10.

These distributed extensions to the Greater Taurus association demonstrate
that an older population is present in the region. The youngest group we find
outside of the well-characterized TMC regions is GT-5 at $\sim$8.5 Myr old,
while all others have ages in excess of 10 Myr. These older regions range from
the subgroups of GT-3, at $\sim$14-19 Myr, to $\mu$ Tau (GT-1), at $\sim$52
Myr old. While their inclusion in Greater Taurus is not necessarily certain
for some of the peripheral subgroups (see Section 5.4), the case of GT-9 is
quite notable as an older, distributed population that is directly connected
to the known elements of the Taurus association. This group consists of two
leaf subgroups, GT-9A and GT-9B, which are contiguous with one another in
spatial coordinates and have overlapping velocity distributions. Subgroup
GT-9A has a near-newborn age solution and contains known Taurus YSOs around
L1517 (e.g., Luhman, 2018), while GT-9B contains a dispersed population with
an older age solution at $\sim$11.5 Myr old, forming a clear link between an
older (age $>$10 Myr), distributed group and the current star formation events
in Taurus. GT-6, GT-4A, and GT-4B also all contain members with known signs of
youth and kinematics previously linked to Taurus (Kraus et al., 2017;
Krolikowski, submitted), and with ages of $\sim$17, 38, and 30 Myr,
respectively, these groups further enforce Taurus’s connection to an older
population.

Our work supports recent studies that have also demonstrated the presence of
this older distributed population, including Krolikowski (submitted), Kraus et
al. (2017). Not all current literature promotes this conclusion, however, such
as Luhman (2018), which concludes that the proposed older members with age
$>$10 Myr are likely not true Taurus members due to differing velocity
distributions. Some of our proposed members have kinematics that make them
potentially subject to this argument, likely requiring a full kinematic
analysis and traceback as discussed in Section 5.4 to achieve a definitive
conclusion on their association with the current star formation in the TMC.
Other old subgroups, however, have clear kinematic overlap with known Taurus
members near the TMC, and those with more discrepant velocities are linked by
previously unknown low-density stellar distributions to the known Taurus
groups, strengthening those connections. The presence of 9 groups and
subgroups in the 10-40 Myr age range, as well as a much older group in $\mu$
Tauri with an age of $\sim$52 Myr, demonstrates that not only is there an
older stellar population surrounding the currently known, essentially newborn
TMC-assiocated groups, but that Greater Taurus may represent the oldest
association in the solar neighborhood that still hosts active star formation.

###  4.2 Greater Orion

Table 5: Clustering within Greater Orion. ORI is the Greater Orion ID, a
number from EOM clustering used for identification. Leaf clusters within an
EOM cluster are given a subcluster ID, which is given in the LEAF column.

ORI | LEAF | Name | N | RA | Dec | l | b |  $D_{sky}$aaOn-sky spatial extent in galactic l/b, in the form of the RMS in major axis $\times$ minor axis, when fit with a bivariate Gaussian. | d | $\mu_{RA}$ | $\mu_{Dec}$ | $V_{T,l}$ | $V_{T,b}$ |  $\sigma_{V_{T}}$bbThe radial extent of the velocity distribution, in the form of the RMS in semi-major axis $\times$ semi-minor axis, when fit with a bivariate Gaussian. | Age
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
|  |  |  | (deg) | (deg) | (deg) | (pc) | (mas/yr) | (km/s) | (km/s) | (Myr)
1 |  |  | 458 | 80.6 | -3.6 | 205.8 | -21.4 | 8.3$\times$3.6 | 316.8$\pm$15.7cchas a member within 1 pc of the search horizon, and therefore members are likely present beyond that limit. | 1.1 | -0.9 | 1.9 | 0.8 | 1.3$\times$1.0 | 13.3$\pm$0.3ddSubstructure present, but mostly coeval.
1 | A |  | 18 | 78.7 | -7.1 | 208.2 | -24.7 | 2.1$\times$0.8 | 288.6$\pm$7.4 | 1.7 | -3.0 | 4.7 | 0.2 | 0.6$\times$0.3 | 14.0$\pm$1.2
1 | B |  | 10 | 80.9 | -1.5 | 203.9 | -20.1 | 1.1$\times$1.0 | 327.7$\pm$1.7cchas a member within 1 pc of the search horizon, and therefore members are likely present beyond that limit. | 1.1 | -0.5 | 1.5 | 1.2 | 0.5$\times$0.3 | 12.3$\pm$1.6
1 | C | 25 Ori | 22 | 81.5 | 1.6 | 201.3 | -18.2 | 2.5$\times$0.4 | 328.4$\pm$2.4cchas a member within 1 pc of the search horizon, and therefore members are likely present beyond that limit. | 1.5 | -0.2 | 1.4 | 1.9 | 0.6$\times$0.3 | 13.2$\pm$0.8
2 |  | Eridanus North$\dagger$$\dagger$footnotemark:  | 82 | 67.8 | -10.2 | 205.9 | -35.7 | 8.8$\times$5.1 | 227.8$\pm$18.2 | 2.3 | -3.8 | 4.8 | 0.4 | 0.9$\times$0.6 | 20.0$\pm$1.8
3 |  | Orion Southeast$\dagger$$\dagger$footnotemark:  | 81 | 87.5 | -8.4 | 213.6 | -17.4 | 4.0$\times$2.6 | 294.1$\pm$8.5 | -0.0 | -0.6 | 0.7 | -0.4 | 1.1$\times$0.6 | 19.0$\pm$0.9

††footnotetext: name is newly assigned by this paper.

Figure 15: HDBSCAN-identified groups and subgroups in Orion, labelled
according to their Orion ID, and shown in RA/Dec, galactic X/Y, and transverse
velocity coordinates. The diamonds and squares correspond to leaf subclusters
of the larger ORI-1 group (yellow circles), while icons are top-level EOM
clusters. Small grey dots represent unclustered Greater Orion members. The
main stars of the Orion constellation are included for reference.

The Orion OB1 association contains an enormous population of stars in terms of
both numbers and extent, while also hosting the Orion Nebula Complex, the
nearest site of active high-mass star formation to the Sun (Hillenbrand &
Carpenter, 2000; Bally, 2008). Substructure has long been recognized in the
region, with Blaauw (1964) dividing the region into four distinct
substructures based on location and photometry, labelled Orion OB1a, b, c, and
d. All of these Orion OB1 subgroups are centered beyond this project’s horizon
at 333 pc, however both Orion OB1a and OB1b have components that are nearer
(Brown et al., 1994; Briceño et al., 2019). The known populations accessible
through this project therefore consist of Orion OB1b, a young ($\sim$5 Myr)
population immediately surrounding Orion’s Belt, and Orion OB1a, a somewhat
older ($\sim$11 Myr) population containing groups to the northwest, such as
the 25 Orionis and HD 35762 clusters (Bally, 2008; Briceño et al., 2019). All
recognized subgroups of Orion OB1a have significant components within our
search horizon, however Orion OB1b has a notably bimodal parallax
distribution, resulting in only the nearest component being visible, a
subgroup given the name Orion OB1b I by Briceño et al. (2019).

Like in Taurus, we also identify a distributed population beyond the known
subgroups in Orion, albeit with much less variety in age. The high-level
region we identify as Greater Orion hosts three large HDBSCAN-identified
groups, with the largest (ORI-1) subdividing into three small subgroups. Their
extents in spatial and velocity coordinates are shown in Figure 15. The mean
transverse velocities in these groups never differ from one another by more
than 4 km s-1, and all groups have tight velocity dispersions of order 1 km
s-1 or below. Particularly on the far edge of ORI-1, these velocity
distributions are found to be nearly identical to those of known young
associations near to and beyond our search radius cutoff, such as 25 Orionis
and the related subgroups Zari-B0 and Zari-B6 defined by Zari et al. (2019).
Some of the more distant members of ORI-1 we identify, including those in
subgroups ORI-1B and ORI-1C, overlap with the near edge of Orion OB1a (Briceño
et al., 2019), with ORI-1C members overlapping directly with known members of
the 25 Orionis and HD 35762 subgroups in spatial coordinates. Other sections
of ORI-1 overlap with the membership of Orion OB1b I, however most of this
group is beyond our search radius limit, leaving the remaining members too
sparsely distributed for HDBSCAN to identify a subgroup. Both ORI-1B and
ORI-1C, along with the rest of the far edge of ORI-1’s main body also appear
cut off by our search horizon, as expected given the distance distributions of
known Orion subgroups (e.g., Briceño et al., 2019). This indicates that our
coverage of known associations in Orion is consistently incomplete.

The lack of overlap between our extended population in Orion and known Orion
subgroups is demonstrated by the DBSCAN-based subclustering analysis of Orion
conducted in Zari et al. (2019), which investigated the full range of
distances that Orion occupies. The only overlap between our members and those
identified with a subgroup by Zari et al. (2019) is at the far edge of our
ORI-1 group (including ORI-1C), which overlap with the membership for Zari-B
in 3-d positions and proper motions. Zari-B6 in particular has a distribution
in position/velocity parameter space that overlaps fully with the membership
of our subgroup ORI-1C, while the near edge of Zari-B0 also overlaps with part
of ORI-1’s far edge (Zari et al., 2019). However, most of Zari-B, as well as
the rest of the groups they identify are well beyond our search radius.

The more nearby populations accessible through this research therefore largely
consist of stars well outside of known groups, residing in lower-density
extensions that methods like the DBSCAN clustering implementation used in Zari
et al. (2019) are not able to recover due to the inflexibility of DBSCAN in
simultaneously handling groups with very different densities and scales. Many
of the stars in these extensions are linked directly to the known Orion OB1a
and OB1b I populations through membership in ORI-1, forming significant
extensions to the region that stretch over 50 pc away from those known
sections in the approximate direction of the sun, and towards Rigel in the
plane of the sky. As we assign them new names: Eridanus North for ORI-2, and
Orion Southeast for ORI-3.

Eridanus North forms a major low-density extension to Greater Orion, spanning
from the near edge of ORI-1 to within $\sim$200 pc of the Sun, with some
members residing more than 150 pc from the known membership in Orion OB1 (Zari
et al., 2019; Briceño et al., 2019). While there is a subtle separation
between the mean transverse velocities in Eridanus North and the majority of
ORI-1 of approximately $\sim$3 km s-1, a small subset of ORI-1 stars at the
near edge have transverse velocities that overlap with those in Eridanus
North. Stars in ORI-1 that have these unique motions form the ORI-1A subgroup,
which has spatial positioning consistent with ORI-1, but transverse motions
more similar to Eridanus North. As such, ORI-1A appears to represent a bridge
between the main Orion Association and the more dispersed extensions into
Eridanus North. The stars in Eridanus North appear notably older compared to
the known populations in Orion, with an age estimated at $\sim$20 Myr.

The other new major group we find in the Orion Complex, Orion Southeast is
identified as UPK 422 in Sim et al. (2019) , but is otherwise undiscussed in
the literature. The group has limited separation from from the main body of
ORI-1 in transverse velocities but much more visible separation in spatial
coordinates. This separation is most evident in sky coordinates, where the
group appears to the Southeast of Orion’s belt, with a significant central
concentration of stars near the border with Monoceros (see Fig. 15). Like in
Eridanus North, we find relatively old ages for stars in this group, with an
age estimate of $\sim$19 Myr.

Despite the large extensions we find to known populations, ORI-1 is found to
have very consistent ages across its three subgroups, all of which are at
$\sim$12-14 Myr and within uncertainties of each other. Consequently, ORI-1
appears to be mainly coeval, and no significant hints of age gradients are
detectable in the region. One minor exception is a small unclustered group of
stars near Orion’s belt associated with Orion OB1b I, where the members have
individual age estimates consistent with the 5 Myr photometric age estimate
for the group from Briceño et al. (2019). Our $\sim$13 Myr age solution for
ORI-1C is consistent with

###  4.3 Perseus

Figure 16: The age distribution and subclustering in Perseus. Ages in RA/Dec
space are shown in the top panel, displaying the clear divide between the
young populations in the west near the Perseus Molecular Could and older
populations to the east. The remaining panels show subclustering in RA/Dec,
galactic X/Y, and galactic transverse velocity. The two subclusters are
labelled in the middle panel, and objects in the older eastern extensions are
marked with icons to distinguish them from the young western populations. The
centres of IC 348 and NGC 1333 are also marked for reference (Ortiz-León et
al., 2018).