Gaia EDR3 Reveals the Substructure and Complicated Star Formation History of the Greater Taurus-Auriga Star-forming Complex

To construct a large, inclusive sample of Taurus, we have compiled a set of 658 objects from multiple sources: an internal list of canonically accepted disk-bearing Taurus members (Kenyon & Hartmann 1995; Luhman et al. 2010; Rebull et al. 2010; Kraus et al. 2011; Kraus & Hillenbrand 2012; Esplin et al. 2014, among others), a compilation of candidate class III members from (Kraus et al. 2017; objects with a membership assessment of "Y" or "Y?"), and recent work on the Taurus census from Luhman (2018), Esplin & Luhman (2019), and Galli et al. (2019). We do not include the new objects from Zhang et al. (2018), as we finalized the sample based on which objects had sources in Gaia DR2. We separate the HBC 358 system included in Galli et al. (2019) into two objects (HBC 358 A and B), as they are separated widely enough to have separate Gaia catalog entries. ¹ We also compile information on binary companions for the canonically accepted Taurus members, adding companions wide enough to have their own Gaia catalog entries, but excluding any companions that are not resolved in Gaia from the presented census. This includes 63 systems that contain unresolved multiples. Binary status, system architecture, and literature references are included in Table 1 for objects where applicable.

For each object in our sample, we queried the Gaia EDR3 catalog (Gaia Collaboration et al. 2021) for sources within 5'' of an object's position. For searches resulting in multiple matches, we used the astrometry and photometry of the sources found to match the object correctly. This process was straightforward for single stars, as nearly all are bright enough to have Gaia measured astrometry, and incorrectly matched sources had wildly discrepant astrometry or photometry. This census crossmatch was first performed with the Gaia DR2 catalog (Gaia Collaboration et al. 2018a), and we used the Gaia DR2 to EDR3 designation-matching catalog as another test to make sure sources were matched correctly.

The crossmatch for Gaia resolved binary systems was more complicated, as the closest match in the Gaia catalog is not necessarily the correct object match. We compared the separation and position angle between the Gaia sources to literature value in order to ensure it is the correct pair of objects. We then assigned a match based on the position angle and Gaia photometry, comparing to literature contrast as an extra check.

In all, 587 objects in our census have Gaia EDR3 source detections (compared to 568 for DR2), with 528 objects having full astrometric solutions (compared to 508 for DR2). The objects with Gaia EDR3 counterparts are listed in Table 1. The 71 objects included in the initial census that do not have Gaia EDR3 sources are listed in Table 2. Figure 1 shows the sky positions of all objects in our census of the greater Taurus region. Of the 104 known binary pairs in our census that have literature separation and contrast values (including pairs within higher-order systems), 37 have source detections for both components, and 26 have full astrometric solutions for both components. All but one of the systems with separations above 1'' are resolved in Gaia EDR3, which we adopt as the resolution limit for Gaia EDR3. Systems with separations smaller than 1'' that are resolved are all close to equal brightness, and there are resolved systems with separations as small as 300 milliarcseconds. We find that a vast majority of wide binary systems that have been identified using various imaging techniques but are unresolved in Gaia have RUWE values above 1.2, which supports the use of RUWE as an indicator of binarity (e.g., Lindegren et al. 2021a).

Zoom In Zoom Out Reset image size

This census is constructed to be expansive, but we err on the side of purity in terms of requiring objects to have some previous confirmation of youth. Our census has inherited the biases and incompleteness of its input catalogs. Many surveys of Taurus have focused near the clouds, introducing a location bias. Infrared and x-ray surveys are more likely to find Taurus members away from the clouds, but the former is biased toward younger, disk-bearing objects, and the latter is biased toward young solar-type stars that are active and bright enough in X-rays to be detected. This results in an incompleteness of the census away from the clouds and at ages older than canonical Taurus, qualifying the conclusions that can be drawn about older, distributed populations of Taurus. Large-scale, consistent searches for young (τ ≲ 50 Myr) stars without any initial phase space preference (e.g., Zari et al. 2018; Kounkel & Covey 2019; McBride et al. 2020; Kerr et al. 2021) are key to unraveling the picture of the greater Taurus region and any as-yet unknown related populations around it.

We have compiled additional data on the objects in our census needed for further analysis presented in this paper. These additional data include spectral types, extinctions, J magnitudes, and radial velocities, which are listed for all objects with those measurements in Tables 1 and 2.

A majority of the spectral types and extinctions were taken from the compilations of Kraus et al. (2017), Luhman (2018), and Esplin & Luhman (2019), although the original reference is adopted where applicable. Spectral type and extinction values determined simultaneously from spectra were prioritized for adoption. Uncertainties on each are adopted from listed references. For spectral types from Esplin et al. (2014), Luhman (2018), and Esplin & Luhman (2019), uncertainties of ±0.25 and ±0.5 subclasses are used for values derived from optical and infrared spectra, respectively (following Luhman et al. 2017, among others). For objects from Slesnick et al. (2006b), a spectral type uncertainty of 0.5 subclasses is assumed from Slesnick et al. (2006a). The uncertainty on A_V from Herczeg & Hillenbrand (2014) is assumed to be 0.25 mag. For A_V derived using photometry from Kraus et al. (2017), we assumed an uncertainty of 0.3 mag. For objects from Luhman et al. (2009a, 2017), Luhman (2018), and Esplin & Luhman (2019), we assume an uncertainty on A_J of 0.1 mag for extinctions derived from infrared spectra following Luhman et al. (2017) and 0.14 mag for extinctions derived from optical spectra and photometry following Luhman (2004) and Luhman et al. (2003), respectively. Luhman et al. (2017) say that the extinction derived from infrared spectra for objects with A_J > 2 are more uncertain, so we adopt an uncertainty of 0.3 mag for those objects. It is unclear if the same methodology is used across these papers, but these are the only available uncertainty estimates. We adopt an uncertainty of 0.45 mag, the mean uncertainty in A_V for objects with reported extinction errors, for objects without literature uncertainty values.

Regardless of the band of the adopted extinction, we have computed extinctions in V, J, and Gaia EDR3 G bands. To convert between V and J, we assume an extinction coefficient of R(J) = 0.72 ± 0.01 (Yuan et al. 2013), which corresponds to A_V /A_J = 4.31 assuming R(V) = 3.1. Error is propagated from both the adopted extinction and R(J). To convert to A_G , we adopt a temperature-dependent relation for A_G /A_V derived using theoretical spectral templates. We determine temperatures for the objects in our census using the adopted spectral types.

Nearly all of our adopted J magnitudes were adopted from the 2MASS point-source catalog (Cutri et al. 2003; Skrutskie et al. 2006). For select objects, we have adopted J magnitudes presented in Esplin & Luhman (2019) from other photometry sources.

We compiled radial velocities from various literature sources for our sample, prioritizing higher-precision literature values over Gaia catalog radial velocities. Fourteen objects have new radial velocities reported derived from observations with the Tull coudé optical echelle spectrograph on the 2.7 m Harlan J. Smith Telescope at McDonald Observatory (Tull et al. 1995). The spectra are reduced with standard procedures implemented in Python, and radial velocities are extracted using broadening functions with theoretical model spectra (Tofflemire 2019). Roughly 70% (415/587) of the objects with Gaia counterparts have radial velocities, with a median uncertainty of 0.4 km s⁻¹. For objects without a literature radial velocity error, we adopt an uncertainty of 1 km s⁻¹.

For each object in our sample with Gaia astrometry, we derived galactic Cartesian XYZ positions by generating 100,000 Gaia positions and parallaxes using the full Gaia covariance matrix. We apply object-specific zero-point corrections to the Gaia parallaxes using the color, brightness, and location-dependent relation from Lindegren et al. (2021b). ² An analysis of eclipsing binaries by Stassun & Torres (2021) finds that these corrections from the Gaia Collaboration are promising, and we include them for consistency with the Gaia team despite their insignificance at the typical Taurus distance. We then convert to galactic XYZ and adopt the median value. This allows us to compute individual XYZ covariance matrices for all objects in our sample, which is important because XYZ is highly covariant along the line-of-sight direction.

We combine the literature compiled radial velocities with Gaia proper motions to compute full 3D kinematics of our sample. We compute galactic UVW velocities for all objects with radial velocities by generating 100,000 Gaia positions and astrometric values using the full Gaia covariance matrix, again adopting the median value and computing individual kinematic covariance matrices. The median radial velocity uncertainty is roughly an order of magnitude larger than the reported projected proper motion uncertainty, which will again introduce a line-of-sight covariance to the UVW velocity data. All Gaia-derived quantities are listed in Table 1.

With this expansive Taurus census, we can robustly identify substructure in the region. We choose to identify clusters in galactic Cartesian coordinates, rather than with astrometric quantities, and excluding any kinematic information. We do not cluster in R.A., decl., and parallax, as they do not form an orthogonal coordinate system for a population of finite extent, and the covariance matrix of a Gaussian with nonzero spread in this space cannot encompass the curvature of the plane of the sky. As many objects do not have radial velocities, and thus do not have full galactic kinematics, we exclude kinematic information from the clustering and will analyze kinematics separately. This should not significantly affect clustering results, especially for well-defined tight spatial groups, and post-clustering kinematic analysis can provide further validation of our group assignments.

To find spatial clusters, we use a Gaussian mixture model (hereafter GMM; McLachlan & Peel 2000; Reynolds 2009), which is a semi-supervised machine-learning clustering model that assumes a data set can be fully described as a combination of multiple Gaussians. Each Gaussian component has its own mean and covariance matrix, allowing for elongated structures to be identified. Group membership is assigned probabilistically, accounting for covariances in the model Gaussians. This is a more accurate and robust model than traditional Euclidean distance-based clustering algorithms, such as k-means clustering (MacQueen 1967; Hartigan 1975), which partitions data into clusters by simply minimizing the square of the Euclidean distances between data points and their assigned cluster means, limiting model group shapes to be spherical in three dimensions. ³

Each component in a GMM has three defining parameters—a mean vector, a covariance matrix, and a weight. The mean and covariances describe the location and shape of each Gaussian component, and the weight describes the relative prior probability for each component, which is essentially the fraction of objects in the fit data set belonging to each component. In a GMM, we can describe the probability distribution of a given data set X as the sum of an ensemble of Gaussians:

$$\begin{eqnarray}&&p(X)=\sum _{i=1}^{K}{w}_{i}\,{ \mathcal N }(X,{{\boldsymbol{\mu }}}_{i},{{\rm{\Sigma }}}_{i}),\end{eqnarray} \tag{ 1 }$$

where K is the number of components in the mixture model, ${ \mathcal N }$ denotes the multivariate normal distribution, μ _i and Σ_i are the mean vector and covariance matrix of the ith Gaussian component, and w_i is the component weight.

With this model, data points are assigned probabilities of being generated from each Gaussian component relative to the other components, rather than labels as in other clustering techniques. In this work, we assign objects the label corresponding to the highest component membership probability. We track the membership probabilities in each group for all objects in our analysis in case there are ambiguous group assignments.

As we are searching in galactic XYZ space, we can only define the GMM with the 528 objects in our census that have full Gaia astrometric solutions. Additionally, we exclude 16 objects that have anomalously discrepant XYZ locations with respect to the extended Taurus region, being too far away (with X > 225 pc); none of these objects would be consistent with the broader Taurus region even accounting for parallax error. While these objects would likely just add additional outlier groups to the GMM fit, they would be extraneous and potentially bias the properties of actual groups related to Taurus on the outskirts of the region. These 16 excluded objects are denoted with an "NM" (for nonmember) in the Group Label column in Table 1.

We make no cuts to the sample based on the quality of an object's astrometric solution. Poor parallaxes could affect how the GMM finds clusters, particularly as the parallax error comprises nearly all of an object's XYZ error and causes the values to be highly covariant in the line-of-sight direction. We address this later quantitatively with detailed fitting of well-defined groups accounting for each object's individual XYZ covariance (see Section 3.2).

To perform the GMM fit to our data set, we use the GaussianMixture module from the open source scikit-learn machine-learning package in Python (Pedregosa et al.2012). We initialize our GMM fit using 100 iterations of k-means clustering. The number of Gaussian components that comprise the model are not fit parameters, but chosen. We calculate and compare the Bayesian information criterion (BIC) for all fits to determine the best fit. The BIC depends on both the output likelihood and complexity of a model, such that the preferred model will have a maximized likelihood while being penalized for the number of Gaussian components. Additionally, the GMM fit is slightly sensitive to the random seed for the k-means clustering initialization.

To ensure we are finding the global BIC minimum, we run fits ranging from 1 to 15 components using 300 different random initial seeds to probe differences with the initialization, resulting in 4500 fits total (300 for each number of components). Figure 2 shows the BIC as a function of the number of GMM components for all fits, along with the minimum and median BIC for each number of components. The global minimum BIC fit has 14 components, which we will adopt as an initial set of groups for further investigation in Section 4.1.

Zoom In Zoom Out Reset image size

One key drawback to the scikit-learn GMM routine is that it does not account for errors on the data themselves. In calculating the galactic XYZ position of the objects in our census, the vast majority of measurement error is introduced by the parallax, which propagates to a strong covariance along the line of sight. For the Taurus region, the X coordinate holds most of this parallax error, with the covariance strongest in the XZ plane. From the 100,000 XYZ samples we calculate for each of the objects in our census, we can compute an individual XYZ covariance matrix for each object, which should be utilized in determining the physical properties of our GMM-found subgroups.

Accounting for the data covariance is important to separate physical elongation from apparent elongation along the line of sight. To incorporate the data covariance information, we perform an MCMC fit of the individual GMM-found Taurus subgroups using emcee (Foreman-Mackey et al. 2013). We model each subgroup as a Gaussian defined by nine parameters: three for the XYZ mean position of the subgroup (${{\boldsymbol{\mu }}}_{\mathrm{mod}}$) and six parameters defining the covariance of the model Gaussian (${{\rm{\Sigma }}}_{\mathrm{mod}}$).

A covariance matrix is a variance matrix that is rotated to introduce correlation between variables. We can decompose a covariance matrix into its eigenvalues and eigenvectors, which are the variances of the (unrotated) variance matrix and rotation matrix, respectively. The six parameters in our model defining the model covariance are three eigenvalues ( λ ) and three angles ( θ _λ ). These parameters can be combined to generate a covariance matrix using the below equation, where ${ \mathcal R }$ is the rotation matrix function and diag is the diagonal matrix function:

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\mathrm{model}}={ \mathcal R }({{\boldsymbol{\theta }}}_{{\boldsymbol{\lambda }}})\,\mathrm{diag}({\boldsymbol{\lambda }})\,{ \mathcal R }{\left({{\boldsymbol{\theta }}}_{\lambda }\right)}^{\top }.\end{eqnarray} \tag{ 2 }$$

With this, we can easily generate random covariance matrices by generating three random angles and three random eigenvalues (greater than zero), without worrying about the need for directly generated random covariance matrices to be positive semi-definite.

We modify the data covariance matrices in two ways to account further for the underreported nature of Gaia astrometric errors. First, we multiply the data covariance matrices by the square of an object's RUWE. This will downweight objects with worse astrometric solutions in the MCMC fit. If an object is an outlier but has large RUWE, it could very well be a bona fide member of the group but would bias the group's parameters from the GMM. Weighting by the RUWE will decrease the effect this outlier would have on the fit parameters. We also include an extra free parameter in our fit, f_err, which is an error inflation term. While including the RUWE acts as an error inflation term specific to each object, f_err accounts for a uniform underreporting of astrometric error, as is known to be the case in Gaia. We do not allow for the case that astrometric errors are overreported, so we require f_err ≥ 1.

It is straightforward to incorporate the modified data covariance matrix into the likelihood calculation because covariance matrices add. Thus, the covariance matrix utilized in the final likelihood calculation for an individual data point, Σ_fit,i , is

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\mathrm{fit},i}={{\rm{\Sigma }}}_{\mathrm{model}}+{f}_{\mathrm{err}}\,{\mathrm{RUWE}}_{\mathrm{data},i}^{2}\,{{\rm{\Sigma }}}_{\mathrm{data},i}.\end{eqnarray} \tag{ 3 }$$

The likelihood function over all N data points is then

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{log}{ \mathcal L } & = & -\displaystyle \frac{1}{2}\sum _{i=1}^{N}{({{\boldsymbol{x}}}_{{\boldsymbol{i}}}-{{\boldsymbol{\mu }}}_{\mathrm{mod}})}^{\top }{{\rm{\Sigma }}}_{\mathrm{fit},i}^{-1}({{\boldsymbol{x}}}_{{\boldsymbol{i}}}-{{\boldsymbol{\mu }}}_{\mathrm{mod}})\\ & & +\mathrm{log}(8{\pi }^{3}\,\det {{\rm{\Sigma }}}_{\mathrm{fit},i}).\end{array}\end{eqnarray} \tag{ 4 }$$

In the MCMC fit, we do not allow the error inflation term to be less than one, nor any of the model covariance matrix eigenvalues to be less than zero, and we restrict the model covariance matrix angles to be between −π and π. We also apply a Jeffreys prior on the covariance matrix eigenvalues:

$$\begin{eqnarray}&&{\mathrm{logPrior}}_{{\boldsymbol{\lambda }}}=-\displaystyle \frac{1}{2}\sum _{n=1}^{3}\mathrm{log}{\lambda }_{n}.\end{eqnarray} \tag{ 5 }$$

All parameters are initialized using the output for each particular subgroup from the GMM fit, and f_err is initialized uniformly between values of one and six. We also run two different types of fits: one in which we fit each of the subgroups individually, and one in which we fit all of the subgroups together. In the latter, we use one global f_err term. This also allows us to compare the fit f_err values for each group individually to a global value, which will inform the degree to which each group is elongated along the line of sight, and any connection between error inflation and further substructure in a group.

Figure 3 shows the results of the MCMC fit on one of our final adopted subgroups, Group C2-L1495, using a group-specific f_err value. The group's spatial extent is shown using the probability density function of the fit Gaussian. The fit from the GMM is more extended than for the MCMC, which is smaller almost exclusively along the line of sight. The MCMC fit is also more centrally concentrated, meaning that the group is more compact than it appears in the GMM. Objects with RUWE above 1.4, which denotes a worse astrometric solution, are highlighted. Interestingly, these points tend to be on the outlying regions of the group along the line of sight. This confirms that the quality of the astrometric solution is significantly affecting the parallax, causing an apparent elongation along the line of sight.

Zoom In Zoom Out Reset image size

The initial conditions and subsequent star formation history of the parent molecular clouds are imprinted in the substructure of the extended Taurus-Auriga region's stellar population. We use the GMM fit described in Section 3.1 to quantitatively identify subgroups in Taurus and assign group membership to our census, which will provide a basis with which to study the spatial, kinematic, and temporal structure of the broader Taurus region. In this section, we analyze the output of the GMM fits, investigate simple properties of the identified subgroups, and finalize the subgroup membership assignments for our census.

Figure 4 shows the XYZ positions of our sample, plotted with colors and markers denoting their membership in groups identified in the initial GMM fit (note that not all 17 groups in the legend are used here, as further group subdivision will be explained below). There appear to be three types of groups identified by the GMM: (1) well-defined groups that are relatively confined in XYZ and mostly reflect the core of the Taurus region (hereafter referred to as core groups and plotted as circles); (2) broad, sparser groups that have some apparent small-scale clustering but are found throughout XYZ space (hereafter referred to as distributed groups and plotted as stars); and (3) groups with very few members (N ≤ 3) that are physically wide (hereafter referred to as outlier groups and plotted as stars).

Zoom In Zoom Out Reset image size

We adopt a naming convention for the subgroups that conveys both the group type (core or distributed) and location. The names of core groups all start with "C," and the names of distributed groups all start with "D," followed by a number for distinction. Each group name has a description or cloud name appended to the end to convey its location within the Taurus region using just text. The outlier groups are named as distributed groups (D5-D7). There have been previous naming schemes for groups in Taurus, such as the Roman numerals to distinguish groups on sky from Gomez et al. (1993) and Luhman et al. (2009a), and the multiple different groupings using Gaia in recent work (Luhman 2018; Galli et al. 2019; Roccatagliata et al. 2020). However, we do not adopt the Roman numeral naming, as the 2D grouping does not map exactly to the 3D grouping. We also do not adopt any names from recent work, as their groups are either too broad, too small, or use a significantly different Taurus census. The definitions of the groups in Taurus are ever-evolving, leading to difficulty in adopting a single consistent naming scheme; we include a cloud name where possible, to best orient the reader to the group locations.

This delineation can be seen in the number density of the groups, which is visualized in Figure 5. This plot shows the number of members in each group as a function of its 1σ ellipsoid volume. There appear to be four separate density regimes: groups with small volume and nearly an order-of-magnitude spread in membership count (Groups C1-L1551, C2-L1495, C5-L1546, C6-L1524, and C9-118Tau), groups with many members and very large volume (Groups D2-L1558 and D4-North), groups in between these two regimes (Groups C3-L1517, C7-L1527-B213, D1-L1544, and D3-South), and groups with very few members across multiple orders of magnitude in size (Groups D5-OutlierCentral, D6-OutlierNorth, and D7-OutlierEast). It is clear that the first of these regimes includes the core groups, the second includes the distributed groups, and the last consists of the outlier groups. For the groups in between, we consider Groups C3-L1517 and C7-L1527-B213 to be core groups, and Groups D1-L1544 and D3-South to be distributed groups. While Group C3-L1517 is extended, it appears to have a central concentration of stars with an extended halo, and Group C7-L1527-B213 appears to have two physically distinct, compact loci. Groups D1-L1544 and D3-South are considered distributed groups, as they are sparser and less well-ordered. These groups with intermediate density may represent the transition from compact group at formation to a dispersing, evolved group; assigning them a core or distributed label is necessary for investigations of further subdivision.

Zoom In Zoom Out Reset image size

The core groups are elongated along the line-of-sight direction (denoted by the background dashed lines in Figure 4), particularly Groups C2-L1495, C3-L1517, C5-L1546, and both visually distinct stellar overdensities in Group C7-L1527-B213. We discuss this apparent elongation in Section 4.2. Groups C1-L1551 and C3-L1517-Halo have halos of objects further from the mean of their respective groups (significantly more so for Group C3-L1517-Halo), but are still relatively compact and do not clearly overlap with any other group. Additionally, Group C7-L1527-B213 appears to be two separate, spatially connected groups with visually distinct loci. The distributed groups, while having some apparent clustering, are mostly spread out across the full XYZ space considered.

To explore whether or not there is further substructure in the region not picked up by the original GMM (such as the apparent subdivision in Group C7-L1527-B213), we ran additional GMM fits on two subsets of our sample, the core and distributed groups, separately. We include the outlier groups with the distributed groups, although there is no meaningful difference if they are included or not. We use the same methodology as with the initial GMM fit.

For the core groups fit, the global minimum BIC is achieved with 11 components. We also look at the minimum BIC 10 component fit, which has virtually the same BIC as the global minimum fit. These two fits are nearly identical, splitting the two visually distinct loci of the original Group C7-L1527-B213 (into the adopted Groups C7-L1527 and C8-B213), distinguishing between the core and halo of the original Group C3-L1517 (split into the adopted Groups C3-L1517-Halo and C4-L1517-Center), and splitting two members off of the original Group C9-118Tau. Group C9-118Tau was already small, but those two objects are distinct in XYZ. Groups C9-118TauE and C10-118TauW include four of the ten members of the 118 Tau association as defined in Gagné et al. (2018). These two groups are apparently associated with the 118 Tau association, and thus have adopted that name. The only difference between the two fits is a subdivision in Group C6-L1524 in the 11 component fit that does not visually stand out. For this reason, and the nearly identical BIC values, we adopt the minimum BIC 10 component fit for the core groups, bringing our total number of core groups identified to 10.

For the distributed groups fit, the global minimum BIC is achieved with six components, although there is significant overlap in BIC values for fits with four to eight components. Inspecting these fits, they preserve the original outlier groups, while adding one or two more small membership outlier groups and slightly rearranging the larger distributed groups. Notably, Group D1-L1544 is kept exactly as it was in the original fit in all inspected subset fits. While there is apparent further substructure in these distributed groups, such as spatial correlations in all groups and overdensities in Groups D1-L1544 and D3-South, the GMM struggles to pick up this clustering due to the broad and sparse nature of these groups. This could be explained by dispersion as clusters age, which would agree with the conclusion that these distributed groups are older than the core Taurus region (Kraus et al. 2017). As such, we adopt the distributed and outlier group properties found with the original GMM.

For our final GMM group membership assignments, we adopt the core groups from the 10 component subset fit, and the distributed and outlier groups from the original 14 component fit. The BIC value of the adopted GMM ensemble is plotted as the horizontal dashed line in Figure 2, and is smaller than the global minimum of the original set of fits, assuring that it is a reasonable fit in terms of membership assignment and complexity. Figure 6 shows the XYZ positions of our sample, along with a legend for the color and marker for each of these adopted Gaussian components (which will be used throughout the rest of the paper). Figure 5 also shows the subdivisions of the core groups found in the adopted GMM ensemble. Various properties of the 17 spatial subgroups are given in Table 3.

Zoom In Zoom Out Reset image size

Table 3. Taurus Group Properties

Group	(α, δ, π) a	Distance	(X, Y, Z)	(U, V, W)	Nmem	Age b	Marker c
	(deg, deg, mas)	pc	pc	km s−1		Myr
C1-L1551	(68.083, 18.0493, 6.909)	144.74	(136.03, 2.09, −49.41)	(15.75, −15.22, −7.38)	44	${1.73}_{-0.2}^{+0.24}$
C2-L1495	(64.3793, 28.1684, 7.72)	129.53	(122.23, 24.27, −35.36)	(16.22, −12.28, −10.72)	55	${1.34}_{-0.18}^{+0.19}$
C3-L1517-Halo	(74.089, 29.5632, 6.356)	157.33	(154.56, 17.98, −23.22)	(16.02, −14.64, −11.09)	36	${2.34}_{-0.38}^{+0.52}$
C4-L1517-Center	(74.1644, 30.3172, 6.444)	155.19	(152.48, 19.25, −21.53)	(14.8, −14.84, −10.47)	26	${2.53}_{-0.5}^{+0.97}$
C5-L1546	(68.7759, 22.8293, 6.22)	160.76	(153.76, 11.66, −45.46)	(16.47, −14.23, −6.9)	32	${2.01}_{-0.29}^{+0.32}$
C6-L1524	(68.0159, 25.2741, 7.788)	128.40	(122.95, 14.53, −34.04)	(15.78, −11.53, −9.49)	72	${1.57}_{-0.17}^{+0.18}$
C7-L1527	(70.2087, 25.7439, 7.051)	141.81	(137.16, 13.92, −33.25)	(15.33, −11.61, −9.5)	22	${2.59}_{-0.78}^{+0.74}$
C8-B213	(66.5258, 26.5957, 6.413)	155.94	(148.54, 22.7, −41.69)	(17.53, −13.08, −6.77)	21	${3.09}_{-0.72}^{+0.92}$
C9-118TauE	(83.8523, 22.845, 9.129)	109.54	(108.83, −7.54, −9.87)	(13.3, −18.94, −8.79)	6	${6.13}_{-1.64}^{+2.63}$
C10-118TauW	(80.4945, 24.8425, 9.485)	105.43	(104.72, −1.08, −12.2)	(13.25, −18.73, −8.71)	3	${41.46}_{-24.36}^{+34.97}$
D1-L1544	(78.3837, 24.1845, 5.949)	168.09	(166.2, −0.12, −25.14)	(19.13, −14.15, −8.81)	24	${3.40}_{-0.72}^{+0.9}$
D2-L1558	(69.7134, 17.3824, 6.784)	147.40	(139.22, −1.7, −48.37)	(15.43, −14.03, −8.64)	55	${3.27}_{-0.36}^{+0.42}$
D3-South	(66.2198, 22.0251, 8.138)	122.89	(115.87, 10.75, −39.49)	(14.01, −6.81, −9.52)	25	${6.22}_{-1.68}^{+1.43}$
D4-North	(66.6245, 25.7332, 7.0)	142.86	(135.99, 19.01, −39.42)	(16.17, −12.66, −9.17)	83	${2.49}_{-0.34}^{+0.35}$
D5-OutlierCentral	(71.5097, 20.0666, 4.724)	211.71	(203.33, 1.53, −58.95)	(16.12, −8.21, −5.14)	3	${17.66}_{-9.61}^{+22.62}$
D6-OutlierNorth	(65.6635, 29.6987, 4.38)	228.30	(217.03, 44.49, −55.12)	(16.82, −7.78, −4.42)	2	${4.71}_{-2.28}^{+8.13}$
D7-OutlierEast	(77.5648, 20.4195, 17.612)	56.78	(55.6, −2.66, −11.2)	(10.04, −5.95, −9.64)	3	${23.23}_{-9.93}^{+34.53}$

Notes.

^a Gaia EDR3, J2016. ^b Uncertainties are from the 16th and 84th percentiles of the age posteriors. ^c The marker and color used for each group in the plots in this paper.

Download table as:  ASCIITypeset image

We recompute membership probabilities for all objects in the sample using these 17 Gaussian components. In the adopted fit, there are 317 members of the 10 core groups, 187 members of the 4 distributed groups, and 8 members of the 3 outlier groups.

Almost all of the core groups have substantial spread in the line-of-sight direction, especially compared to their spread in the plane of the sky, with Groups C2-L1495, C4-L1517-Center, C5-L1546, and C8-B213 being the most elongated. Star formation often produces filamentary structures (Hartmann 2002; Molinari et al. 2010), at least until those structures relax into clusters or disperse into the field (Baumgardt & Kroupa 2007; de Grijs & Goodwin 2008). As such, it is not unexpected that we would identify groups that are elongated along a particular axis. However, it would be highly coincidental for these groups to be truly elongated near exactly along the line-of-sight direction. To determine whether or not this line-of-sight elongation is real, we fit the core groups while accounting for the covariance in the astrometric data and the underreporting of error in the Gaia EDR3 catalog. We do this by multiplying the data covariance matrices by both the square of their individual RUWE values, to capture the higher uncertainties for sources with poor fits, and by an error inflation term that is the same for all objects in a particular fit. We performed two types of MCMC fits: one on each group individually with their own error inflation term, and one on all of the core groups simultaneously with a global error inflation term.

The global fit finds a median ${f}_{\mathrm{err}}={1.12}_{-0.08}^{+0.11}$, which means that Gaia errors are underreported by around 12%. This agrees with the roughly 10% error underreporting found by Gaia Collaboration et al. (2018b), which is roughly the same as in EDR3 (Fabricius et al. 2021). El-Badry et al. (2021) looked at wide binaries in EDR3 and used their reported parallaxes and corresponding uncertainties to assess the degree to which the errors are underreported. They provide a relation to compute the expected parallax error inflation term for a single star given its G magnitude, which for the objects in the core groups produces a median ${f}_{\mathrm{err}}={1.06}_{-0.03}^{+0.08}$, in agreement with, but smaller than, the results of our global fit. The higher error inflation we find is likely attributable to additional noise introduced by the presence of unresolved binary companions, through orbital motion or photocenter jitter as the two stars vary in brightness independently over time.

Figure 7 shows the f_err posterior distributions from the MCMC group fits, including fits with a global error term and individual group-wise error terms. When the groups are fit individually, all groups except Group C6-L1524 find an f_err larger than the global fit, ranging from 1.2 to 1.9. The f_err values for Groups C4-L1517-Center, C5-L1546, and C8-B213 do not agree with the global f_err. This is not surprising, as these three groups are all highly elongated along the line of sight, confirming that the parallax error underreporting is linked to this elongation. Group C9-118TauE favors a much larger f_err (2.4) with a long tail toward even larger values because it only has six members, and thus the fit is ill-constrained and results in a large Gaussian fit because it has little structural knowledge to inform the fit. While we do not include them in the global fit, we did perform individual fits on the distributed groups. We exclude them because they are significantly broader than the core groups and do not have obvious elongation along the line of sight. Groups D2-L1558, D3-South, and D4-North all have f_err smaller than the global fit, while Group D1-L1544 finds a larger f_err, due to it having a cluster of outlier points.

Zoom In Zoom Out Reset image size

While the global fit provides insight into the catalog-wide error inflation present in Gaia EDR3, we adopt the individual fits for the purposes of analyzing the structure and size of groups. The individual fits are better at accounting for outliers in a group, in particular outliers that have relatively good RUWE values, as an increased f_err will be able to downweight those data in the fits when RUWE cannot. We interpret the individual f_err values not solely as direct measurements of the error inflation, but as a convolution of the error inflation term with properties of the individual groups. This includes further substructure, the frequency of undetected binary systems, and the spatial distribution of group members with poorly constrained astrometric solutions.

To examine the extent to which our subgroups are affected by error inflation along the line of sight, we compute the standard deviation for each group both along the line of sight and in the plane of the sky. To do this, we sample from the Gaussian function defining each group from both its GMM and MCMC fit. We then project the XYZ positions of these samples onto the line-of-sight direction at the group's mean location. The norm of that projected vector is an object's line-of-sight coordinate, and the line-of-sight variance is then the variance of all group members' line-of-sight coordinates. We compute the plane of sky vector as the difference between the XYZ position and line-of-sight vector, and again the plane-of-sky coordinate is the norm of this vector. This plane-of-sky projection is a one-dimensional displacement for a coordinate in a two-dimensional plane. To get the plane-of-sky variance, we find the displacement that encloses 39.3% of the samples, which corresponds to the 1σ distance for a 2D Gaussian.

Figure 8 shows the line-of-sight and plane-of-sky standard deviations for both the GMM and MCMC fits of the core groups. The line-of-sight variance decreases for all groups, confirming that Gaia error elongated the groups along the line of sight. This is reinforced by the fact that the plane-of-sky variance is virtually unchanged in all groups. All groups also fall closer to the one-to-one line with the MCMC fit results compared to the GMM fit results, meaning they are closer to equal aspect ratio (although they may have structure in the plane of the sky). The distributed groups do not have substantial change in either their line-of-sight or plane-of-sky variances, confirming that their large spatial extents are physical and not apparent.

Zoom In Zoom Out Reset image size

To validate the shrinking along the line-of-sight, we additionally ran fits arbitrarily setting f_err = 1 to ensure the inflation term is not solely responsible for shrinking the spatial extent of these groups. The fits with a free f_err are slightly smaller, particularly along the line of sight, which is expected since the additional uncertainty should allow for a smaller Gaussian to encompass a group's members. Otherwise, the fit Gaussians are largely similar in shape and orientation, and this assures that the bulk of the line-of-sight shrinking is due to the distribution of the lowest-uncertainty objects being least elongated, as quantified through the RUWE-corrected data covariances.

As we are accounting for the parallax effect, we can be confident any remaining evidence of structure is real and not apparent. Groups C2-L1495 and C8-B213 are still larger in the line of sight than the plane of the sky, by roughly a factor of two. Group C8-B213 does not have many members, which could be affecting the fit quality, but Group C2-L1495 is very well-populated. This indicates that Group C2-L1495 and Group C8-B213 are in reality elongated along the line of sight. These groups are associated with the clouds L1495 and B213, which have been suggested to have clumpy or layered structure in the line-of-sight velocity projection by Hacar et al. (2016). Our work adds that the cloud is likely spatially structured in the line of sight as well, and that Groups C2-L1495 and C8-B213 are the result of star formation along filaments in the cloud that have collapsed earlier than the present gas.

If the groups we have identified in Taurus are from distinct formation events, we expect them to have small internal velocity dispersion, with larger differences between the groups. The relative velocities of the groups can also shed light on the star formation history of the region, such as whether or not some subsets of the groups formed together. Figure 9 shows the bulk galactic UVW velocities for all subgroups, which are presented in Table 3, along with the standard deviation in their mean velocity vectors. We only include objects that have RV errors less than 0.5 km s⁻¹, and exclude objects that have compiled RVs from LAMOST, as they are subject to systematic uncertainties larger than our criterion (Guo et al. 2015; Luo et al. 2015; Zhong et al. 2015; Kiman et al. 2019; Luo et al. 2019; Zhong et al. 2019). Two of the outlier groups (D6-OutlierNorth and D7-OutlierEast) each have only one member with a good RV, while the third outlier group only has two. Therefore, we do not analyze their kinematics further.

Zoom In Zoom Out Reset image size

The core groups, except for Groups C9-118TauE and C10-118TauW, which will be discussed later, have bulk motions within 4 km s⁻¹ of each other in all three kinematic dimensions. There is coherence present, linking Groups C1-L1551, C5-L1546, and C8-B213 together, and Groups C2-L1495, C3-L1517-Halo, C4-L1517-Center, C6-L1524, and C7-L1527 together, although the latter set of groups has more dispersion. While there is some distinction in V, nearly all of this separation occurs in W. This parallels the spatial clustering: Groups C1-L1551, C5-L1546, and C8-B213 are all more distant and lower in Z than the other set of groups. Groups C6-L1524 and C7-L1527 are particularly close in both kinematics and spatial position. The kinematics of Groups C3-L1517-Halo and C4-L1517-Center agree, which is consistent with their spatial coincidence. They also have a dearth of precise RVs compared to other groups, leading to a larger uncertainty in their bulk motions, which is most readily apparent in their large extent in U. Group C3-L1517-Halo also has a larger velocity dispersion than Group C4-L1517-Center, which is consistent with being a dispersed halo of objects around Group C4-L1517-Center. The RV uncertainties dominate the errors on the bulk kinematics of these groups. The U velocity carries most of the RV information, and for all core groups has the largest error of the three velocity components.

The two core groups that do not follow this categorization are Groups C9-118TauE and C10-118TauW, which are significantly removed from where the main Taurus groups sit in velocity space, particularly in V. This is not surprising, as these two groups are significantly removed from the main Taurus group spatially, are similar to each other spatially, and older. Their kinematics, however, are not discrepant with being linked to Taurus. On average, they are separated from the other core groups by roughly 50 pc in space and 6.3 km s⁻¹ in velocity. This is in line with expectations from Larson's Law (∼5–7 km s⁻¹, depending on the formulation), and is similar in scale to spatial and kinematic substructure seen in the Sco-Cen association (Rizzuto et al. 2011). It is possible that this small group was formed in a slightly earlier epoch of star formation linked to the recent and ongoing star formation in the main Taurus region. If Groups C9-118TauE and C10-118TauW are a part of the greater Taurus ecosystem, ongoing and future searches for young stars in the region between them and the core groups should have velocities that fall in between their bulk motions.

The distributed groups are less kinematically homogeneous than the core groups, and do not all reside in the same general location in velocity space as the core groups. The distributed groups have larger typical velocity errors than the core groups. This is particularly evident in their larger V and W error bars, which are much less affected by the RV error, which dominates the uncertainty budget. The larger velocity errors are indicative of a larger velocity dispersion, which can be explained if these groups are dispersed collections of coeval, potentially older stars. Groups D2-L1558 and D4-North are both kinematically close to the core groups, sitting in roughly the center of the core groups' UVW distribution. This bolsters the idea that these two groups comprise stars either born in an earlier epoch of star formation near the Taurus region that have dispersed or that have formed in isolation concurrently with the core groups. In fact, while Groups D2-L1558 and D4-North are spatially distinct from each other (especially in Y), they can visually be joined together to create one large group of objects interspersed throughout Taurus with roughly the same kinematics. This distinction in Y could be due to these groups originating from two main sites of ongoing star formation: near Group C1-L1551 and clouds L1551 and L1558 for Group D2-L1558, and in the main canonically defined Taurus complex near Groups C2-L1495, C8-B213, and C6-L1524, and clouds L1495, B213, L1521, and L1524 for Group D4-North. This is analogous to Groups C1-L1551, C5-L1546, and C8-B213, which are spread across tens of parsecs in the Y direction, yet are kinematically coherent.

Groups D1-L1544 and D3-South are more distinct, with Group D1-L1544 removed from the other groups in U and Group D3-South removed in V. Neither of these groups is particularly spatially distinct, although Group D1-L1544 is separated from the other distributed groups. Group D3-South directly overlaps Groups D2-L1558 and D4-North. Even without kinematics, the GMM identified the snake-like structure of Group D3-South crossing through Groups D2-L1558 and D4-North. This structure is further highlighted by the fact that Group D3-South's plane of sky spatial standard deviation is nearly three times larger than that along the line-of-sight. Coupled with its distinct kinematics, Group D3-South is a distinct subgroup of Taurus with a filamentary or planar structure unlike nearly all of the other identified groups.

Stars form in molecular clouds that are turbulent, which is important in the star formation process as it leads to and regulates fragmentation, helps set the initial mass function, and is the seed of a stellar association's velocity field. Larson's Law is an empirical relation that relates the size of molecular clouds to their velocity dispersion, which is set by the turbulent power spectrum (Larson 1981). The velocity structure of our sample can directly test if Larson's Law applies to the groups we have found in Taurus in describing the turbulent nature of the region, and how much evolution occurs in the velocity structure after stars form and decouple from the gas. To do this, we compute the velocity structure function (VSF; McKee & Ostriker 2007) of our sample, calculating the pairwise velocity difference between objects as a function of their 3D spatial separation. We calculate VSFs for multiple subsets of our data, such as using only objects that are members of core groups or distributed groups together.

Computing VSFs using just the positions and velocities of the objects themselves does not account for the uncertainties on their 6D phase space coordinate, largely from parallax and RV error. There are also covariances introduced, as different VSF points can have the same start or end point, so uncertainty in one position can introduce a covariance on the velocity structure information of N − 1 baselines across the N-body system. To account for the uncertainty and covariance introduced in the VSF computation, we generate 1000 samples of galactic position and kinematics for each object included in a VSF. We draw these realizations of position and velocity from their Gaia astrometric solution and covariance matrix, and our literature-compiled RV and corresponding error. To avoid inflation of the VSF by observational uncertainties, we exclude all objects that do not meet the RV quality condition described in Section 5.1. With this, we produce 1000 different instances of the VSF encapsulating the underlying uncertainty in the velocity structure of the Taurus stellar population.

Figure 10 shows three sets of VSFs for the core and distributed groups separately: one in which all objects in each type of group (core or distributed) are compared to each other ("all point comparison"), one in which objects are only compared to other objects in their specific subgroup ("intra-group comparison"), and one in which objects are only compared to other core or distributed objects not in their specific subgroup ("inter-group comparison"). There are differences between each of the subsets of the sample that are highlighted by the fits shown, particularly between the core and distributed group subsets.

Zoom In Zoom Out Reset image size

The VSF data follow power laws, as expected from Larson's Law and formation in a turbulent molecular cloud, and we fit their slope and intercept to quantify the turbulent structure of the region. For a particular VSF, we fit each of the 1000 VSF instances separately, and then adopt values and uncertainties from the resulting slope and amplitude parameter distributions. We exclude all VSF points with 3D spatial separations below 1 parsec, as that is the typical distance uncertainty, and any points below that will be effectively smeared out and hold no precise structure information.

For the fit, we first bin each VSF into equally spaced logarithmic bins in 3D spatial separation. The bin edges are computed using the astropy implementation of Knuth binning, which is a Bayesian estimation of the optimal binning that is robust to the underlying data distribution (Knuth 2006). We bin using the corresponding data value-only VSFs (the adopted phase space coordinates of each included object, no sampling); we use the same bins on all VSF instances. Within each bin of 3D spatial separation, we calculate the median velocity difference, and fit the resulting binned VSF with a power-law function. We exclude any bins with fewer than five data points in them, to reduce biasing the fit from outlier points with small number statistics, particularly at the smallest and largest scale separations. We adopt the median of the 1000 fits as the VSF parameter values and the 16th–84th percentile range as the 1σ parameter uncertainties. For some of the fits, as will be described below, we fit multiple power laws over different spatial separation regimes.

The distributed group VSFs are accurately described with single, well-behaved power laws for each of the three cases. The slope for the intra-group VSF is larger than the slope for the all-point VSF, which in turn is larger than the inter-group VSF. Both the intra- and inter-group VSF slopes, however, agree with the all-point VSF within 1σ. As the distributed groups are large and sparse, there is not a clear spatial scale delineation between the inter- and intra-group comparisons. The intra-group VSF power-law slope is similar to the original formulation of Larson's Law (σ_v = 1.1 L^0.38; Larson 1981). Assuming that the stars have retained memory of the parent molecular cloud's turbulent structure, this would indicate that the groups formed from a cloud with subsonic turbulence. The inter-group VSF has a shallower slope, which could result from the different distributed groups forming from different clumps in the larger molecular cloud, such that they have retained the velocity structure within their clump but did not inherit any velocity structure beyond that level.

The core group VSFs are significantly more complex. The all-point VSF is a clear broken power law, with a change in amplitude and slope at roughly 17 pc. The smaller-scale regime is dominated by the intra-group VSF, and the larger-scale regime is dominated by the inter-group VSF. To determine if there is further underlying structure, or if there is simply a change of power law at that scale, we must inspect the separate inter- and intra-group VSFs.

The intra-core-group VSF is a steady, single power law across its entire separation range, and has a slope nearly identical to that of the small-scale regime of the all-point VSF. The inter-core-group VSF is quite different from the intra-core-group VSF, with two clear regimes: a nearly flat VSF below 17 pc and a power law above 17 pc with a slope nearly identical to that of the all-point large-scale VSF. The inter-core-group VSF has a large jump in velocity dispersion, suggesting a significant energy injection into the entire core group region at a ∼17 pc scale. There are few data points in the intra-core-group VSF above that scale, so it is hard to conclude if the energy injection is evident there as well. The inter-core-group core VSF is essentially flat at 3 km s⁻¹ between 4 and 17 parsecs, which means that even where individual core groups are closer to each other, there is no velocity structure between their members.

We also interpret the difference between the inter- and intra-core-group VSFs at small scales to further validate the core group assignments. If these group identifications were spurious, we would expect the two VSFs to be identical at small scales. However, we see that the velocity difference between objects is smaller when they are members of the same group than when they are not at a given separation below 10 pc. This is another confirmation that group membership is meaningful with respect to possible association with the other groups. There are many fewer objects at the smallest scales in the inter-distributed-group VSF, which makes it hard to determine if this is true for the distributed groups as well, although we expect that to be the case (and do see some apparent flattening at the smallest scales).

There are two major ways in which the VSF of a collection of stars can evolve over time, one that steepens the relation and one that makes it shallower. Objects at smaller separations, while having lower velocity dispersion, will move fractionally further away from each other than objects that are at large separations. This will result in the objects at small scales "catching up" to the objects at large scales, steepening the slope of the VSF. On the other hand, groups of stars can dynamically relax for as long as they remain bound, losing memory of the initial velocity structure of the region and producing a shallower slope (asymptotically approaching a flat VSF) where there is no pairwise spatial and kinematic correlation.

This process occurs over the relaxation timescale, which is given as:

$$\begin{eqnarray}&&{t}_{\mathrm{relax}}=\displaystyle \frac{N}{6\mathrm{log}\tfrac{N}{2}}\,\displaystyle \frac{R}{{\sigma }_{v}},\end{eqnarray} \tag{ 6 }$$

where N is the number of stars in the group, R is the radius of the group (taken as the half-mass radius), and σ_v is the velocity dispersion; the second fractional term in the equation is the crossing time. For our groups, we can only estimate a relaxation timescale because their membership may not be complete and our number count is a lower limit. Also, relaxation may begin before molecular cloud material is expelled, leading to a deeper potential well with faster typical stellar velocities. We assume that our groups are mass-mixed, and that the half-mass radius is equivalent to the half-count radius. To calculate the half-count radius, we take 10,000 samples using each group's mean and covariance matrix output from the MCMC spatial fits using an individual f_err term. We then determine the 1D distance at which half of the samples are enclosed. We take the velocity dispersion to be the standard deviation in 1D velocity separation for objects in a group that pass the RV quality cut from 5.1.

We can interpret the VSFs presented here in the context of these two processes in order to draw conclusions about the nature of the groups. The slope of the intra-core-group VSF is much shallower than that of the intra-distributed-group VSF. The core groups are much more compact than the distributed groups, with calculated relaxation times ranging from 0.5 to 3 Myr for the core groups and 4.3 to 9 Myr for the distributed groups. With these short relaxation times, combined with the typically assumed Taurus age of less than a few million years, we conclude that the individual core groups have begun to dynamically relax while sufficient molecular cloud material remained to keep them bound, resulting in a shallower VSF slope. If the distributed groups are young and formed as they are seen now—sparsely and across a large spatial extent—they have not had enough time to relax, and will more closely retain the velocity structure of the parent molecular cloud. If they are older, evolved, dispersing equivalents of the core groups, they may have already dynamically relaxed to some extent, before gas expulsion unbound the groups and their dispersal led to a steeper VSF slope.

It is harder to apply these processes to the inter-group VSFs, in particular for the core groups, as they are distributed through the region in spatial clusters, and thus have no well-defined crossing time, unlike the distributed groups, which are regularly spaced throughout the region. The inter-core-group VSF's steeper slope is in agreement with the intra-distributed-group VSF's slope, implying that the original velocity structure of the cloud has not been lost on the larger scale between groups. It is unclear why the inter-distributed-group VSF has a shallower slope, but it perhaps is an indication that the distributed groups are less related to each other than the core groups are to themselves. If they formed at different times, unlike the core groups, we do not expect them to behave as though they formed from the same molecular cloud.

The slope of the velocity structure function of stars in Taurus is shallower than the currently accepted typical Larson's Law slope of 0.5 from Solomon et al. (1987). Qian et al. (2018) measured the velocity structure function of molecular cores in Taurus identified in Qian et al. (2012) over projected separation scales from 1 to 10 pc. They find that the molecular core VSF follows a power law with a slope between 0.5 and 1/3, with a steeper slope below 3 pc and a shallower slope above 5 pc. This is in agreement with the roughly 1/3 slope we find for the VSFs we believe trace the initial spectrum of the stellar population, although we find much higher velocity dispersions at all separations than in the molecular core VSF. Stars inherit the velocity structure of their parent molecular cores, except exhibiting higher-velocity dispersion, perhaps from early dynamical interactions when the stars decoupled from their parent cloud. The agreement between the VSFs of these two different types of objects implies that the turbulent power spectrum in Taurus may be shallower than other molecular clouds.

To derive ages for the Taurus subgroups, we construct a Bayesian framework for calculating ages from isochrones for a collection of stars. To calculate the likelihood of a star (or system) falling at a specific HR diagram or magnitude–magnitude diagram position, we compute the underlying probability density function of a synthesized stellar population of a given age. We draw each synthetic star's magnitude from an isochrone and broaden the ensemble's magnitude–magnitude sequence using a kernel density estimator (KDE) to generate a smooth distribution. With the KDE likelihood, we compute a posterior distribution of age. We finally explore the posterior space using a Markov chain Monte Carlo with emcee (Foreman-Mackey et al. 2013) to find the best-fit age and its uncertainty.

For this work, we use Gaia EDR3 G and 2MASS J magnitudes. We do not use Gaia G_BP or G_RP magnitudes, as they both can be particularly sensitive to the presence of binaries, and they have less flux and fewer scans than G. Due to the common presence of disks, we choose 2MASS J over K, as the disk excess will add additional K-band flux above the photospheric level, which we do see bias our fits toward younger ages if K magnitudes are used.

We first generate the ensemble of model stars of the given age with masses drawn from a three-power Kroupa IMF, as presented in Equation (2) of Kroupa (2002). We then interpolate masses and age on an isochrone grid to get G and J magnitudes for each model star. Here, we use a combination of the BHAC15 and MIST isochrone grids (Baraffe et al. 2015; Dotter 2016; Choi et al. 2016). The MIST grid is only computed down to 0.1 M_⊙, while the BHAC15 grid is computed down to 0.01 M_⊙. We adopt the MIST grid above 0.1 M_⊙ and the BHAC15 grid below. There is a slight discontinuity at the transition between the two models (0.1 M_⊙), but it is smaller than the typical measurement uncertainty. The MIST grid has been updated to include the Gaia EDR3 bandpasses, while the BHAC15 grid hasn't, but the differences between the MIST DR2 and EDR3 magnitudes are much smaller than typical uncertainties. We therefore have chosen to use the BHAC15 grid's DR2 colors in interpreting the observed EDR3 photometry.

We have also implemented a binary population synthesis to account for unresolved multiples. Following Sullivan & Kraus (2021), we assign a binarity flag to each model star using a mass-dependent multiplicity fraction shown in Table 4. We then draw a separation for each binary system using the mass-dependent log normal distributions described in Sullivan & Kraus (2021). This prescription covers objects with masses between 0.1 and 2 M_⊙; we adopt the 0.1 M_⊙ distribution for objects with masses below this range and the 2 M_⊙ distribution for objects more massive than this range. We then draw a mass ratio from a power-law distribution whose index is mass-dependent. We use the same piecewise linear interpolation described in Sullivan & Kraus (2021) to generate these mass-dependent power-law distributions, with one exception: we account for the dichotomy seen in the mass ratio distribution at high primary masses. De Rosa et al. (2014) find that A-star binaries with separations smaller than 125 au have a more uniform mass ratio distribution than wider-separation systems. We therefore replace the highest mass point in the piecewise interpolation with a power-law index value of −0.5 for systems with separations below 125 au. For the binary systems, we compute magnitudes for each component from the interpolated isochrone grid.

Depending on the binary system's separation, we then combine the individual object magnitudes to get an unresolved system magnitude. We assume a common distance of 145 pc to convert all physical separations into angular separations, and use different angular resolutions for each band. For Gaia G, we consider all systems with separations below 1'' unresolved, while we treat wider systems as fully resolved with the secondaries having their own model Gaia sources and include them separately in the model population. For 2MASS J, we consider all systems with separations below 25 to be unresolved. For systems with separations between 25 and 4'', we treat the system as resolved in 2MASS but do not include the secondary as its own source separately in the model population (Cutri et al. 2003). For systems with separations larger than 4'', we include the secondary as a separate object in the model population. Synthetic stars that only have one magnitude, exclusively the G magnitude of secondaries in binaries with separations from 1'' to 4'', are excluded from the ensemble.

To create a continuous likelihood distribution from this synthetic stellar population, we convolve the (M_J , M_G ) synthetic population data with a KDE to calculate a smooth likelihood function in magnitude–magnitude space. We use a Gaussian kernel with a bandwidth of the typical uncertainty in the de-extincted magnitudes for our data, which is roughly 0.2 mag (taking into account magnitude, distance, and extinction uncertainties). We can then evaluate the output of the KDE at a specific combination of absolute G and J magnitude to compute a likelihood.

Within the MCMC fit itself, we use a linear-flat prior on age between the bounds of our isochrone grid, excluding ages younger than 0.5 Myr and older than 100 Myr. While 0.5 Myr is the lower limit set by the isochrone grid, it also represents the expected duration of the Class 0 + 1 stages, and hence is the point when stars become easily observable (Evans et al. 2009). For MCMC sampled ages within the grid and the ensemble of objects, the value of the posterior distribution is simply the combination of the likelihood for each of the individual objects:

$$\begin{eqnarray}&&\mathrm{log}P(\mathrm{age}| G,K)=\sum _{i=1}^{N}\mathrm{log}{ \mathcal L }({G}_{i},{K}_{i}),\end{eqnarray} \tag{ 7 }$$

where N is the number of stars in the population, ${ \mathcal L }$ is the smoothed likelihood function from the KDE, and the log-prior term is excluded because it evaluates to zero. To speed up computation, we precompute a grid of KDEs at intervals of 0.01 Myr and use the KDE with an age closest to the sample age within the MCMC in the likelihood calculation. While this discretizes the likelihood space, it is much smaller than the typical uncertainty in the fit age, and the difference between KDEs at this age spacing is minimal. This framework is flexible and can be extended to include other parameters for future investigation, such as the IMF itself, binary demographics, or an age spread within a subgroup.

We exclude from all fits objects that: (1) have A_V > 5, (2) have been previously identified as Class 0 or I, and (3) have spectral types earlier than F5. These objects have less reliable absolute photometry, and are often outliers that negatively affect the quality of the fit.

We have applied the Bayesian formalism of Section 6.1 to each of the groups identified in Section 4.1, determining their best-fit ages and corresponding uncertainties. We have also fit all core groups together and all distributed groups together to assess the ages of the region at large.

Figure 11 shows the fit output for one of our identified groups: Group C5-L1546, which is a core group of 29 objects near the central areas of ongoing star formation. The data are shown in both magnitude–magnitude and color–magnitude space along with the KDE of the best-fit median age of ${2.01}_{-0.29}^{+0.32}$ Myr. The posterior distribution on age, and spatial distribution of the group are also shown. The KDE encompasses the data points well, showing that it accurately accounts for the measurement error and captures the spread in the data. The good match between the group's magnitude–magnitude spread and the KDE broadened solely by observational errors implies that there is not a need for an age spread. This is not true of all groups, however, and will be discussed later. We also find the upper uncertainty to be larger than the lower uncertainty for all fits, although it is more pronounced at older ages. This correctly reflects that the magnitude–age gradient becomes shallower toward older ages, eventually resulting in a tail toward older ages in the posterior beyond a few Myr. While there are substantial uncertainties in the absolute ages, even as much as factor of two (e.g., Naylor 2009; Pecaut et al. 2012; Malo et al. 2014; Feiden 2016; Rizzuto et al. 2020), relative ages are measurable to fairly high precision, barring any serious systematic issues shared by many members of a group. A precision relative age measurement is particularly possible with sufficiently large group membership, which is achieved even in the lesser-populated groups such as Group C5-L1546.

Zoom In Zoom Out Reset image size

Figure 12 shows the posterior age distributions from the Bayesian formalism for all groups. The adopted ages and uncertainties for each group are presented in Table 3. In general, the distributed groups are older than the core groups, except for Group C9-118TauE; this is unsurprising, as Group C9-118TauE is tightly defined like a core group but separated from the cloud complex. Groups that have broad, vaguely multimodal posteriors, such as Groups C4-L1517-Center, C7-L1527, and C8-B213, suffer from having fewer objects than other groups. We do not show the outlier groups (Groups C10-118TauW, D5-OutlierCentral, D6-OutlierNorth, D7-OutlierEast), as they all have three or fewer members, which results in poorly constrained ages. There are only modest inter-group age spreads for each type of group (core and distributed), except for Group D3-South, which is significantly older than the other distributed groups, like Group C9-118TauE is in relation to the other core groups.

Zoom In Zoom Out Reset image size

Groups C3-L1517-Halo and C4-L1517-Center, which are spatially coincident and appear to be a central core and dispersed halo around it respectively, have ages that are in good agreement. We ran a fit combining the members of these two groups, which forms a well-defined sequence in the CMD. This results in an age that falls between the reported ages for the two groups, with a well-constrained posterior distribution. Along with their kinematic similarities, we conclude that these two groups are linked to the same local star-forming event. Future surveys of these combined subgroups, with a more complete subcensus, could investigate population properties, such as the mass distribution, to search for radial dependencies in low-mass star-forming events.

The age uncertainty is strongly dependent on the number of objects included in the fit. However, we find that the distributed groups have slightly larger uncertainties than the core groups even when accounting for the number of objects included in a fit. We conclude that the increased uncertainty is a result of the distributed groups having age spreads within them. This is consistent with the idea that the distributed groups are conglomerations of two populations of stars: one from previous epochs of star formation that have dispersed, and one from recent and isolated star formation throughout the region.

To further quantify the age difference between the core and distributed populations, we performed "super-group" fits, which use the Bayesian formalism separately on all members of core groups together and all members of distributed groups together. The posterior distributions of those fits are shown as the first two distributions in Figure 12. We find a core group population age of 1.80 ± 0.08 Myr and a distributed group population age of ${3.17}_{-0.21}^{+0.23}$ Myr. We emphasize that these ages are systematically uncertain in the absolute sense. However, in the relative sense, the distributed group population is nearly twice as old as the core group population, with an age 1.5 Myr older, which is well outside of uncertainties. We find that the distributed group population has a larger age uncertainty than the core group population. This is partially due to the fact that there are fewer objects in the distributed than the core "super-group" fits (163 versus 266). Even after accounting for the number of objects in the fit (assuming uncertainty goes as $\sqrt{N}$), however, the distributed group age uncertainty is twice that of the core group population. We interpret this to mean that the distributed groups have more of an age spread than the core groups.

The age spread is evident in the distributed population's CMD, which is quite broad, and shown in the upper left panel of Figure 13. In the CMD, there are two apparent populations with different ages, and these populations are demarcated by spectral type. The M-dwarf population, and a handful of earlier objects, agree well with the fit age of the whole population; this is not surprising, as the age fit is driven by objects closer to the peak of the IMF. There is a tight sequence of objects with spectral type earlier than M that is older than the rest of the population. The earlier spectral type preference of this older subpopulation is likely due to selection effects: they have been discovered through surveys that preferentially find earlier-type objects. To investigate this older subpopulation, we performed a separate age fit on objects in the distributed group population that have spectral types earlier than M and are below the 3.17 Myr isochrone by more than 1σ in M_G .

Zoom In Zoom Out Reset image size

The age posterior of this older subpopulation is shown in the upper right panel of Figure 13. We find a best-fit age of ${13.34}_{-2.22}^{+3.62}$ Myr, with the corresponding isochrone plotted in the CMD in Figure 13. The posterior is multimodal, which is likely due to higher sensitivity to outlier objects since there are only 36 data points in the fit. Still, the median age visually matches the CMD sequence exactly. The bottom two panels in Figure 13 show the spatial and kinematic distributions of this older population. These objects are preferentially closer than the rest of the Taurus groups, which matches the finding from Kraus et al. (2017) that the distributed population is a collection of older, foreground objects. These objects have a wide spread in kinematics, but have an overdensity at V velocities larger than the rest of the Taurus region.

The spatial and kinematic properties of this older population closely match those of Group D3-South, which largely comprises these older objects. Group D2-L1558 also contributes multiple members to this older population, but it is a relatively insignificant fraction of its full membership, and half of its older objects are not located near the kinematic locus of Group D3-South. Additionally, the spatial distribution and bulk motion of Group D2-L1558 is relatively unchanged if you exclude these older objects. We conclude that Group D2-L1558 is a robust and large group of objects of varying age distributed throughout the Taurus region. The other distributed groups contribute negligibly to this older population.

The core groups are newly formed, relatively tight clusters of stars that are still near their places of birth. The distributed groups are collections of objects of varying ages, from newly formed to over 10 Myr, that are related to and distributed throughout the Taurus region. They are a combination of older objects that formed in clusters that have since dispersed, younger objects that have been ejected from their birth site, and in the largest fraction, objects of mixed age that simply formed in isolation from the areas of intense and clustered star formation. The ensemble of core groups could evolve to look like the distributed groups: a conglomeration of stars with similar spatial and kinematic properties that formed concurrently but have since dispersed and erased the small-scale substructure with which they were born.

Finally, we investigate any trends between the ages of the Taurus subgroups and their galactic positions and velocities. To do this, we fit linear trends to the core groups' ages and bulk position or velocity points. Figure 14 shows the median group ages as a function of their mean XYZ galactic positions. We do not include the outlier groups C10-118TauW, D5-OutlierCentral, D6-OutlierNorth, or D7-OutlierEast in the plot. There are no significant nonzero trends in Y or Z. The X trend, however, is highly significant at 3σ above no slope. This means that, in the canonical Taurus region, age increases with distance. This is not consistent with the idea that the star formation in the canonical Taurus region was triggered from the direction of the foreground older population, but actually from the other direction. While they are not included in the fit, the other distributed groups (D1-L1544, D2-L1558, and D4-North) do generally follow this trend. We find no significant age trends with UVW velocities.

Zoom In Zoom Out Reset image size

As an outlier in age, location, and kinematics, observed again in this work, the older foreground population has not been considered to have relation to the main Taurus region. However, we conclude that none of those things preclude a connection between the older population and the ongoing star formation in the canonical Taurus clouds. While it is kinematically distinct from the core Taurus region, it is only slightly more so than would be expected from Larson's Law, within velocity errors. It is also entirely plausible that the population was birthed from a distinct molecular cloud with its own bulk motion, but was a part of a larger network of gas in the galactic neighborhood from which the canonical Taurus clouds formed. The cloud from which this population formed could have become dense enough to trigger star formation before the canonical Taurus region, and as a likely similar low-mass star-forming event, it never produced a supernova to trigger star formation in the surrounding areas. The members of Group D3-South that are in the older foreground are slightly closer in X than the rest of Group D3-South, showing that there is at least some gradient toward younger ages between the older foreground population and the canonical Taurus region.

Instead, the canonical Taurus region became dense enough to start forming stars from a trigger further away in X, independent of the older foreground. If it were a supernova, it could be the source of the energy injection into the core group population at a scale of ∼15 pc seen in the inter-core-group VSF in Figure 10.