SPICY: The Spitzer/IRAC Candidate YSO Catalog for the Inner Galactic Midplane

The YSO selection is largely based on IRAC photometry from GLIMPSE (Benjamin et al. 2003; Churchwell et al. 2009) and related surveys that used similar observing strategies and data reduction methodologies. ⁸ These include GLIMPSE I (31,184,509 sources), GLIMPSE II (19,067,533 sources), and GLIMPSE 3D (20,403,915 sources), the Vela-Carina (2,001,032 sources) (Majewski et al. 2007; Zasowski et al. 2009), Cygnus X (3,913,559 sources) (Beerer et al. 2010), and SMOG (2,512,099 sources) (Winston et al. 2019) surveys. Spitzer's observations of the Galactic Center (Stolovy et al. 2006) were included in the GLIMPSE II Catalog. We use only the Spitzer photometry obtained during the cryogenic mission, which includes four mid-IR bands centered at 3.6, 4.5, 5.8, and 8.0 μm. We omit the GLIMPSE 360 data from the warm Spitzer mission that includes only the 3.6 and 4.5 μm bands.

The GLIMPSE team designed their reduction pipeline to provide reliable point-spread function (PSF) fitting photometry in crowded fields with spatially varying nebular emission (Benjamin et al. 2003; Kobulnicky et al. 2013)—conditions that are common in IRAC images of star-forming regions. They provide two source lists for each survey, the "Catalog," which is more reliable, and the "Archive," which is more complete. ⁹ Following, Povich et al. (2013), we use the "Catalog" photometry. We make no additional cuts on quality flags, but certain flags are discussed in Appendix A. The Spitzer/IRAC images have PSFs with full widths at half maximum of 166 at 3.6 μm, 172 at 4.5 μm, 188 at 5.8 μm, and 198 at 8.0 μm. This is significantly better than the ∼6'' resolution provided by the WISE survey over a similar wavelength range (Wright et al. 2010), giving IRAC a distinct advantage in crowded fields in the Galactic midplane. The catalogs from the GLIMPSE team also include many stars that are missing from the Spitzer Enhanced Imaging Products (SEIP) due to GLIMPSE's better treatment of mid-IR nebulosity. PSF photometry is also more accurate than SEIP aperture photometry in regions with variable backgrounds (Fang et al. 2020).

The GLIMPSE I, II, 3D, Galactic Center, and Vela-Carina survey observations consisted of (2–5) × 1.2 s integrations at each positions, while the Cygnus X and SMOG surveys used a 0.4 + 10.4 s high dynamic range mode. Given that Cygnus X and SMOG are deeper than the rest of the data, we impose uniformity by omitting sources in these fields that are either brighter or fainter than the limits for the main GLIMPSE survey. ¹⁰ The magnitudes of our selected YSO candidates range from [3.6] = 8.3–14.9 mag, [4.5] = 7.3–13.7 mag, [5.6] = 6.4–12.9 mag, and [8.0] = 5.5–12.2 mag (1%–99% quantiles), and the median photometric uncertainties are 0.062, 0.073, 0.075, and 0.060 mag in these bands, respectively.

Even in the full GLIMPSE catalogs, the presence of red sources and several catalog artifacts can be seen in color–color diagrams (Figure 1). The highest concentration of sources have colors close to 0 (expected for normal stars without IR excess), but numerous sources form a distribution extending to the upper right in these plots, which is made up of both YSOs and other red mid-IR sources (e.g., evolved stars, galaxies, etc.; see Appendix B). In the [4.5]–[5.8] versus [5.8]–[8.0] diagram, a streak can be seen extending from the origin to the upper left. This streak appears to be related to erroneous photometry in the 5.8 μm band for a low fraction of the GLIMPSE sources, and it extends from the origin because this is where the source density is highest. Another prominent feature in this diagram is a finger-like structure extending upward at [5.8]–[8.0] ≈ 1.6, which we attribute to polycyclic aromatic hydrocarbon (PAH) emission (Section 5.7).

Zoom In Zoom Out Reset image size

Near-infrared JHK_s photometry from the Two Micron All Sky Survey (2MASS; Skrutskie et al. 2006) is already included in the GLIMPSE (and extensions) data products. 2MASS has a spatial resolution of ∼2'', which is comparable to the Spitzer/IRAC PSF. For our sample, 2MASS is nearly complete down to J ∼ 15.4 mag, H ∼ 14.2 mag, and K_s ∼ 13.0 mag, with median photometric uncertainties of 0.038, 0.040, and 0.035, respectively. While these limiting magnitudes correspond well with the limits of the GLIMPSE surveys, in practice YSOs are often found in regions of high interstellar reddening, where 2MASS may not be deep enough to detect NIR counterparts of red GLIMPSE sources.

Deeper NIR catalogs with higher spatial resolution are available from the United Kingdom Infra-Red Telescope (UKIRT) Infrared Deep Sky Survey (UKIDSS; Lawrence et al. 2007) and the Visible and Infrared Survey Telescope for Astronomy (VISTA) Variables in the Vía Láctea survey (VVV; Minniti et al. 2010) for the northern and southern portions of the Galactic plane, respectively, with overlap around the Galactic center. These catalogs are deeper than 2MASS, but are saturated for brighter sources. We use the UKIDSS catalog from the Galactic Plane Survey (Lucas et al. 2008) and the averaged VVV photometry for multiple epochs from the VVV Infrared Astrometric Catalog (VIRAC DR1; Smith et al. 2018). For both deeper NIR surveys, the photometry extends to J ∼ 19 mag, H ∼ 18 mag, and K_s ∼ 16 mag, with formal photometric uncertainties <0.01 mag.

IRAC and UKIDSS/VVV were cross-matched using a 1'' match radius in TOPCAT (Taylor 2005). Photometric measurements from UKIDSS/VVV were omitted if they did not have the flag mergedClass = −1 (stellar) or if they had magnitudes J < 11, H < 12, K < 10.5 (UKIDSS) or J < 12.5, H < 13, K_s < 11.5 (VVV), for which saturation effects start to affect photometry. The higher spatial resolutions of the NIR catalogs mean that it is possible for multiple NIR sources to be associated with individual IRAC sources; however, examination of IRAC+UKIDSS matching by Morales & Robitaille (2017) revealed that the NIR flux is usually dominated by a single counterpart.

In our analysis, we perform the YSO candidate selection independently on cross-matches of IRAC+2MASS, IRAC+UKIDSS, and IRAC+VVV, and the results of these separate selections are merged.

The Gaia mission (Gaia Collaboration et al. 2016), in its second data release (Gaia DR2; Gaia Collaboration et al. 2018), has provided optical broadband photometry (Evans et al. 2018) for the whole sky along with exquisite astrometric measurements (Lindegren et al. 2018) for more than 1.3 billion stars. These data on their own can be used for selecting possible pre–main-sequence stars (e.g., Zari et al. 2018), but for our study we use them as ancillary data in order to better understand the parallax (ϖ) and proper motion (μ_α*, μ_δ ), distributions of the IR-excess selected YSO candidates. From the YSOs candidate list (Section 4), 33% have Gaia counterparts with the full five-parameter astrometric solution (Lindegren et al. 2018). A match rate below 50% is expected because many YSOs are enshrouded by dust and thus not optically visible.

Longer-wavelength photometry is available from both Spitzer's MIPS Galactic Plane Survey (MIPSGAL; Carey et al. 2009; Gutermuth & Heyer 2015) at 24 μm and the WISE All-Sky Data Release (Wright et al. 2010) at 22 μm. In the Galactic plane, AllWISE is affected by high numbers of spurious sources, particularly in the longer-wavelength bands, so we follow the catalog cleaning recommendations from Koenig & Leisawitz (2014) and apply the signal-to-noise and χ² quality cuts on the profile-fit photometry given in their Equations (1)–(4). For MIPSGAL, we use the "Catalog" instead of the "Archive." The WISE photometry, and to a lesser extent the MIPS photometry, is strongly affected by crowding and nebulosity in the regions around the Galactic plane that we are investigating, leading to fewer reliable source detections in these areas. The application of the quality cuts from Koenig & Leisawitz (2014) leave visible holes in the spatial distribution of WISE sources surrounding the clusters of YSO candidates that we identify with the IRAC photometry. Only 30% of our IRAC candidates have counterparts in the 24 μm MIPS band, and 8% have reliable 22 μm WISE photometry. Because of the unavailability of longer-wavelength data for the majority of our sample, we do not use these bands for identifying candidates, but only for post-selection examination of the sample. For both catalogs, we used a cross-match radius of 12.

YSOs identified in earlier studies of star-forming regions within GLIMPSE (and extensions) can be used to train a classifier to find similar types of objects. We use YSOs identified as part of the Massive Young Star-Forming Complex Study in Infrared and X-ray (MYStIX; Feigelson et al. 2013), in addition to a similar, earlier study of the Carina Nebula (Townsley et al. 2011). From the combined lists, we have included probable YSOs from the Carina Nebula, NGC 6611, M17, NGC 6530, M20, NGC 6357, NGC 6334, RCW 38, RCW 36, and DR 21, ranging from ∼0.7 to 2.7 kpc in distance. Povich et al. (2011, 2013) performed the IR excess detection for these projects using Spitzer data, based on a strategy that included SED fitting of both reddened stellar atmospheres and the Robitaille YSO models, color cuts to remove certain types of contaminants, spatial filtering to remove objects that are not clustered, and visual examination of SEDs.

We select all objects from the MYStIX IR-excess catalogs classified as both a YSO (Cl = 0) and a probable member (Mm = 1). More recently, Gaia DR2 has become available, so for the subset of sources with Gaia parallaxes (∼30% of the sample), we refine the sample further by removing any source with a parallax that is discrepant from the median parallax of the group by >2 times the reported parallax error.

In addition to the aforementioned studies of multiple regions and large areas, GLIMPSE has been used in hundreds of papers about individual (or several) star-forming regions. A few representative examples of these include Zavagno et al. (2006), Watson et al. (2008), Povich et al. (2009), Dewangan & Ojha (2013), Samal et al. (2014), Mallick et al. (2015), and Povich et al. (2016).

In the Galactic midplane, many background stars are affected by high levels of foreground extinction, so any source that is either insufficiently red or whose red colors can plausibly be explained by reddening alone is dropped from further scrutiny.

Cuts on IRAC colors and color uncertainties can remove many objects that have no chance of being selected as reliable IR excess objects. We apply the rules recommended by Povich et al. (2011, 2013), decreasing the number of sources in our sample by a factor of ∼10. All retained sources must be detected in at least four out of the seven IR bands, two of which must be 3.6 and 4.5 μm. Sources are kept if there is the suggestion of IR excess in the [3.6]–[4.5] color using the criterion

$$\begin{eqnarray}&&[3.6]-[4.5]-0.408\gt \mathrm{error}([3.6]-[4.5]),\end{eqnarray} \tag{ 1 }$$

where "error" denotes uncertainty in color, calculated by adding the photometric uncertainties for the two bands in quadrature. The value 0.408 is the expected reddening of this color with A_V ≈ 30 mag of extinction. Sources are also kept if they have have photometric measurements in the 5.8 and 8.0 μm bands and either meet both the criteria

$$\begin{eqnarray}&&\left|[4.5]-[5.8]\right|\gt \mathrm{error}([4.5]-[5.8])\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&[5.8]-[8.0]\gt \mathrm{error}([5.8]-[8.0]),\end{eqnarray} \tag{ 3 }$$

or

$$\begin{eqnarray}&&\left|[4.5]-[5.8]\right|\leqslant \mathrm{error}([4.5]-[5.8])\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\left|[5.8]-[8.0]\right|\leqslant \mathrm{error}([5.8]-[8.0]).\end{eqnarray} \tag{ 5 }$$

These rules, optimized from experience with GLIMPSE data, ensure that determination of IR excess is based on more than just the 8.0 μm band, which can occasionally give a spuriously bright measurement.

We fit the JHK+IRAC SEDs of the remaining sources with reddened Castelli & Kurucz (2003) stellar atmosphere models, using the Indebetouw et al. (2005) extinction law. The fitting procedure takes into account the statistical photometric uncertainties on the data, which happen to be of similar size for both 2MASS and IRAC photometry. The UKIDSS and VVV data sets provide JHK photometry for many objects that were not detected in 2MASS, allowing many more sources to be included. However, the statistical measurement uncertainties for most UKIDSS and VVV sources are far more precise than for 2MASS or IRAC. Given that we are mostly interested in detecting deviations from a reddened stellar atmosphere model in the IRAC bands and we want similar selection performance for each data set, we rescale all UKIDSS and VVV error bars that are smaller than the median 2MASS error bars to be equal to the median 2MASS error bars. The sources that are poorly fit by the reddened stellar photosphere, with χ² per data point >4, comprise the target set for our random forest classifier.

Overall, these pruning steps leave 319,251 2MASS+IRAC, 188,701 UKIDSS+IRAC, and 257,334 VVV+IRAC sources with possible IR excess as inputs to our classification step below.

The training data includes both MYStIX IR-excess sources (Section 2.4) that we label "YSO" and sources unlikely to be YSOs that we label "contaminant" (see discussion of contaminants below). Although lists of members are more complete in some of the nearest star-forming regions (e.g., Ophiuchus or Taurus; see Evans et al. (2009b) and Luhman (2018)), we choose to use MYStIX because these massive star-forming complexes may better represent the regions we expect to probe in the Galactic midplane, i.e., at greater distances, with higher extinction, and in more extreme environments. Furthermore, many of the MYStIX regions lie within the survey region of GLIMPSE and its extensions, meaning that homogeneous data products are available for both the training and target sets.

Contaminants can include both sources that occur in star-forming regions (e.g., nonstellar sources such as nebular knots and shocked emission) and sources that are smoothly distributed on the sky (e.g., AGB stars and galaxies; see Robitaille et al. (2008) and Gutermuth et al. (2009)). Refer to Appendix B for discussion of these objects. Povich et al. (2011, 2013) used a variety of techniques to remove these objects from their catalogs, including SED fitting of Robitaille et al. (2007) models, color cuts, and visual inspection. We apply the term "field object" to any object within the field of view of these studies that was not classified as a probable young star by either IR or X-ray criteria.

To enlarge our sample of contaminants, we identify several fields near the Galactic plane that have no signs of star formation (Appendix B.3) and label objects in these fields as non-YSOs. These fields were selected to include lines of sight at multiple Galactic longitudes, in the midplane and up to several degrees above or below it, and with different amounts of Galactic extinction.

Training sets are generated separately for each combination of NIR+IRAC data due to the differences in NIR filters. For 2MASS+IRAC, the training set contains 2865 YSOs, 3436 field objects in the MYStIX fields, and 7718 other field objects for the 2MASS+IRAC data set. For UKIDSS+IRAC these numbers are 919, 2000, and 1128, and for VVV+IRAC they are 1266, 1459, and 2595, respectively.

The distributions of training set objects in color space are discussed in Appendix C. The full IRAC catalog of ∼5 × 10⁷ sources includes a few outliers in region of color space that are not well-sampled by either the labeled YSOs or labeled non-YSOs in the training set. (The limits used to identify these outliers are given in Appendix C.) Given that we have little basis to assign such objects to either category, we opt to be cautious and do not include these objects in our final YSO list.

When data are combined across multiple catalogs, it is almost certain that missing data will occur, as is the case here. Figure 2 depicts the missing pattern for the 2MASS+IRAC, UKIDSS+IRAC, and VVV+IRAC training sets, from which only 57%, 46%, and 29% of objects, respectively, have complete information. About 20% of the 2MASS+IRAC objects are missing three colors at once, and JHK_s are often missing together. The VVV+IRAC data set has the most missing data, with 57% of the objects missing at least three bands. While UKIDSS+IRAC is the most complete, more than 35% of their rows have at least two missing colors. Thus, a naive removal of rows presenting missing values would throw away a non-negligible amount of valuable information.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (3 panels)

Tarball of figure set images

As a final pre-processing step before training our YSO classifier, we employed a multiple-copula imputation. This decomposes joint probability distributions into their marginal distributions and a function, the copula, that couples them (Nelsen 2010). Copulas have been used previously in astronomy—for example, to construct likelihood functions for weak lensing analysis (Sato et al. 2011; Lin et al. 2016) and to infer bivariate luminosity and mass functions (Andreani et al. 2018). Previous tests suggest that this method outperforms other popular approaches, such as multiple imputation via chained equations (van Buuren & Groothuis-Oudshoorn 2011) and Amelia (Honaker et al. 2011), in terms of bias and coverage, especially in cases where the variables are not normally distributed (Hoff 2007). The underlying idea of copula imputation is to derive conditional density functions of the missing variables given the observed ones through the corresponding conditional copulas, and then impute missing values by drawing observations from them. Finally, the choice of performing imputation before training the random forest models has been previously assessed by other studies (Jaeger et al. 2020), which have shown it to reduce the variance in model error estimate, without any detectable change in precision. The imputation method was implemented using the sbgcop package (Hoff 2018) within the r language (Core Team 2019). Copulas were fit simultaneously to both training and target data sets.

Overall, the imputed data preserves the coverage of the original data set (Appendix C). Nevertheless, we do not advocate for the use of these colors in other contexts. They are treated here as nuisance parameters to enable classification of the entire data set.

Decision trees are learning algorithms that resemble the natural flow of human decision making. At each node of the tree, the algorithm randomly selects one feature, and based on the distribution of training data, determines the decision boundaries that best separate different classes. Training objects are then propagated along their branch of the tree to the next node, where a new feature is selected. The process is repeated until the tree reaches a predetermined depth or until all objects in a leaf belong to the same class; for a more detailed description, see Rokach & Maimon (2014). This basic concept has given rise to successful algorithms in many different fields. However, a single decision tree trained on an entire data set is prone to overfitting, presenting low-accuracy results whenever faced with data not used in training. This problem can be overcome by randomizing different stages of the tree construction and combining many independent estimators in a more robust classifier. This type of approach belongs to the wider class of ensemble models.

Ensemble methods (e.g., Sagi & Rokach 2018) are regression algorithms, constructed from the combination of many weak classifiers that, when considered together, provide a more robust estimate than any of their constituents. Random forests (Ho 1995; Breiman 2001) are one such algorithm, composed of many decision trees, each constructed independently. The final classification is determined via majority vote, considering all trees in the forest. In this context, the probability of being a YSO is approximated by the percentage of trees in the ensemble voting for a YSO candidate—we call this probability estimate the "YSO score." Random forests have been successfully used to classify YSOs in smaller-scale studies, including with missing data imputation (e.g., Ducourant et al. 2017; Melton 2020).

We construct the YSO random forest classifiers using the following covariates: J − H, H − K_s , K_s –[3.6], [3.6]–[4.5], [4.5]–[5.8], [5.8]–[8.0], and the 3.6 μm band magnitude. To mitigate the influence of class imbalance upon tree optimization, we employed an upsampling method during the training process. Overall, this does not have a perceptible effect on the classification performance, but it does push the decision boundary probability between classes closer to 0.5. Each model was independently trained for 2MASS + IRAC, UKIDSS + IRAC, and VVV + IRAC data sets. The random forest was employed using the caret R package, with 1500 trees, which was sufficient to guarantee stable solutions. As a sanity check, we tested few other regression models (including generalized additive models, support vector machine, gradient boosting machines, and conditional random forest), but no significant difference in the final YSO candidate set was found. This suggests that random forests (or other typical nonlinear classifiers) capture the data complexity well enough without the need for highly complex models.

We assessed the performance of the classifiers by using validation tests in which we partitioned the labeled objects into a training set (80%) and a test set (20%).

Figure 3 (far left) shows true-positive rate (TPR) as a function of false-positive rate (FPR), a curve often called the receiver operating characteristic (ROC). The area under the curve (AUC) can range from 0 to 1, where random guessing gives an area of 0.5, and higher values indicate better performance. Our classifiers all have AUC > 0.9. A threshold on the YSO score of p = 0.5 gives high TPRs and low FPRs for each classifier, meaning that this threshold is effective for distinguishing between the YSO and non-YSO classes. Other performance measures are listed on the plot, including accuracy (ACC), i.e., the number of true positives plus true negatives divided by the total population, and precision (PRC), i.e., the number of true positives divided by the number of true positives plus false positives. The confusion matrix (Figure 3, center left) shows few misclassifications.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (30 panels)

Tarball of figure set images

Variable importance (Figure 3, center right) is evaluated via out-of-bag samples, which consist of random samplings of the data that are left out of each tree. These are calculated by measuring variations in the prediction error when the out-of-bag data are permuted solely among a specific color, leaving the others unchanged. The process is then repeated across all trees. The final result is a measure of the incremental error for a given color when compared with the unperturbed colors for all the 1500 trees over the entire forest.

Figure 3 (far right) displays a partial dependence plot (PDP; Greenwell 2017) for [3.6]–[4.5] color (plots for other variables are included in the figure set). PDPs are useful for visualizing the relationship between individual features and the response while accounting for the average effect of the other predictors in the model. The shape and steepness of the curves are indicators of the predictor's relative influence. Note the sharp behavior of [3.6]–[4.5], one of the best indicators of YSO candidates.

On the J versus J − H diagram (Figure 9, left), both candidate YSOs and probable contaminants occupy a triangular region of color–magnitude space, where the upper edge of the triangle is approximately parallel to the reddening vector. The YSO candidates are densest around J ∼ 15.5 mag and J − H ∼ 1.3 mag, whereas the probable contaminant distribution is multimodal, with one peak just blueward of the peak of the YSO candidates, and another strip of stars along the upper right edge of the triangle. The stars in this strip, which are more luminous than the typical YSO candidate with the same J − H color, lie in the region of the diagram that would be occupied by reddened post–main-sequence stars.

On the [4.5] versus [4.5]–[8.0] diagram (Figure 9, right), the YSO candidates form a smooth distribution ranging from the bright limit at [4.5] = 6.5 mag to ∼14 mag, where sensitivity declines, with the peak of the distribution at [4.5] ∼ 12.2 mag and [4.5]–[8.0] ∼ 1.2. A ∼ 1 Myr old (pre–)main-sequence star with a mass in the range 0.4–10 M_⊙ at a distance of ∼1–2 kpc would have an unreddened photospheric magnitude 9 ≲ [4.5] ≲ 14 mag (Bressan et al. 2012)—approximately where were find the bulk of the YSO candidates. The probable contaminant distribution peaks at both bright and faint magnitudes. The bright contaminants form a band that tends to be bluer than the YSOs in [4.5]–[8.0], while the faint contaminants tend to have redder [4.5]–[8.0] colors.

The gray lines on the [4.5] versus [4.5]–[8.0] diagram were defined by Gutermuth et al. (2009) to separate dusty AGNs from YSOs in their studies of nearby star-forming regions. Although many of the faintest 4.5 μm sources in our sample have been classified as probable contaminants, the region defined by Gutermuth et al. (2009) for selecting AGNs does not appear to separate our classes well. This apparent discrepancy may arise because Gutermuth et al. (2009) examined deeper Spitzer surveys of relatively nearby star-forming clouds at higher Galactic latitudes, where more AGN are expected to be detected, whereas GLIMPSE is less sensitive to this type of contaminant. Furthermore, GLIMPSE includes more distant star-forming regions in which legitimate YSOs will present fainter observed [4.5] magnitude distributions.

Figure 10 shows the distributions of sources in J − H, H − K_s, and H–[4.5]. On the JHK_s diagram, we include a representative isochrone for ∼1 Myr unreddened stellar models. Most of the objects are shifted to the upper right from this curve, in the approximate direction of the reddening vector. However, the distribution of the YSO candidates spreads to redder H − K_s colors, which would be expected for stars with K_s -band excess. Objects with very red J–H > 5 colors are largely classified as contaminants.

On the H − K_s versus K_s –[4.5] diagram, we show both a 1 Myr isochrone for (pre–)main-sequence stars and a 1 Gyr isochrone that also includes post–main-sequence stars. The red-giant branch extends upward to stars with redder H − K_s colors than the (pre–)main sequence, allowing these groups to be better separated. On this plot, the reddening vector points to the upper right. If we consider a line parallel to the reddening vector, starting from the tip of the asymptotic giant branch (as shown by the gray line in the figure), we would expect that many of the stars lying above this line could be evolved stellar contaminants. This is consistent with what the classifier finds; most objects above this line are classified as probable contaminants, while the candidate YSOs are more abundant below this line. The slope of the Indebetouw et al. (2005) reddening vector is not precisely parallel to the upper edge of the source distribution—this may arise due to systematic uncertainties in the reddening law, or it could be a property of IR colors of highly obscured evolved stars.

Figure 11 shows four projections of sources in IRAC color–color space. On the [3.6]–[4.5] versus [4.5]–[5.8] diagram, the YSO candidates are smoothly distributed, with a peak in density around [3.6]–[4.5] ∼ 0.5 and [4.5]–[5.8] ∼ 0.4, and a tail that extends up and to the right. In contrast, the contaminant distribution peaks slightly bluer in [3.6]–[4.5], and the distribution appears bifurcated, with some sources being redder in [3.6]–[4.5] while others are redder in [4.5]–[5.8]. This bifurcation may be related to the types of contaminants. For example, the contaminants to the upper left roughly correspond to a region of color space identified by Gutermuth et al. (2009) (and indicated by the gray boundary) as containing sources produced by knots of shocked emission from H₂, while the contaminants to the lower right were associated with knots of PAH emission. Both areas outlined by Gutermuth et al. are dominated by objects that we classify as probable contaminants, but the edges of the YSO distribution also overlap these boundaries.

On the [4.5]–[5.8] versus [5.8]–[8.0] diagram, the peak density of YSO candidates is redder in both colors than the peak density of probable contaminants. The objects with the most extreme [5.8]–[8.0] colors are nearly all classified as contaminants. These lie in a region of the diagram identified by Gutermuth et al. (2009) (gray boundary lines) as being dominated by unresolved star-forming galaxies. This diagram also includes a finger comprised of both YSO candidates and contaminants, extending to high [4.5]–[5.8] values ranging from ∼2 to ∼3.5, but with [5.8]–[8.0] colors in a restricted range (1.5 ≲ [5.8]–[8.0] ≲ 2.6. Previously published YSO catalogs (e.g., Gutermuth et al. 2009; Rebull et al. 2011) have included a few YSOs in this region of color space; however, the high number of sources identified when examining the entire inner Galactic midplane makes this feature appear much more pronounced. These stars have colors similar to the PAH nebulosity found in star-forming regions (Povich et al. 2013); however, visual inspection of a sample of these sources suggests that the majority are bona fide point sources in all four IRAC bands.

In the bottom two panels of Figure 11, the peaks of the YSO candidate distributions are redder than the peaks of the probable contaminant distributions for each IRAC color. Nevertheless, while the objects with reddest [3.6]–[4.5] tend to be YSO candidates, the objects with most extreme red [4.5]–[8.0] or [5.8]–[8.0] colors are almost all classified as contaminants.

Figure 12 shows the J − H colors, which are the most sensitive to extinction, versus the IRAC colors, which are the most sensitive to IR excess. In J − H, the peaks of the density distributions are slightly redder for the contaminants than for the YSO candidates, but in IRAC colors, the peaks are significantly redder for the YSO candidates than the contaminants. In both cases, the objects with most extreme red colors tend to be classified as contaminants. However, some of the reddest IRAC sources do not appear on these plots because they lack J-band magnitudes.

On all these color–color plots, the blue ends of the distributions are artificially truncated by the selection rules imposed to ensure that the IR excesses are real. Thus, our catalogs will not be sensitive to certain classes of YSOs, including some YSOs with anemic disks or pre–main-sequence stars without disks.

Less than half the YSO candidates are optically visible: for example, Gaia DR2 detects ∼36,000 of them, which comprise ∼30% of the entire sample. The candidates detected by Gaia tend not to be as red in the mid-IR as other candidates (e.g., [3.6]–[4.5] ≲ 1 and [4.5]–[8.0] ≲ 1.2).

Figure 13 shows a Gaia color–magnitude diagram for the visible YSO candidates. Absolute G-band magnitudes, computed using Gaia parallaxes ϖ, are plotted against Gaia G − RP colors. Only sources with signal-to-noise ϖ/σ_ϖ > 3 are included, meaning that the sample of 7686 sources is small in comparison with the total number of YSO candidates. Nevertheless, this sample is useful for evaluating whether the optically bright candidates have properties consistent with pre–main-sequence stars.

Zoom In Zoom Out Reset image size

On the Gaia color–magnitude diagram, we show isochrones for young stars at several ages, ranging from 1 to 50 Myr, from the Bressan et al. (2012) models. We also indicate the effects of reddening, which would shift points down and to the right on this diagram. The wide Gaia bands mean that the effect of reddening depends on the spectrum of the object, so we show three approximate reddening vectors for three colors; more discussion of how this affects selection of pre–main-sequence stars can be found in Herczeg et al. (2019) and Kuhn et al. (2020). Nearly all the candidates lie above the 50 Myr isochrone, and the majority also lie above the 1 Myr isochrone, which is consistent with most of these candidates being very young pre–main-sequence stars.

When photometry is available in the MIPS 24 μm band (or the W4 band at 22 μm), it can be useful for corroborating classifications based on IRAC. For example, the SED at ∼24 μm tends to be more steeply declining for AGB stars, where IR excess is produced in hot dusty winds, in contrast with YSOs' relatively cooler disks and envelopes.

Figure 14 shows the candidate YSOs and contaminants in J − K_s versus [4.5]–[24] colors. These colors may be useful for distinguishing between AGB stars and YSOs, because the typical AGB star has a precipitous rise in the JHK bands followed by a drop in the mid-IR. The figure shows that the YSO candidates tend to be in the middle of the [4.5]–[24] distribution. The objects with [4.5]–[24] ≲ 2.4 are almost all classified as probable contaminants; however, there is a red tail to the probable contaminant distribution, with a high percentage of objects redder than [4.5]–[24] ≳ 7 also being considered contaminants. An examination of the spatial distribution (not shown) of the probable contaminants in this red tail reveals that many of them are nonclustered, higher-latitude objects in the Galactic bulge. Of the contaminants with low [4.5]–[24], the distribution of J − K_s ranges from ∼0 to ∼9, extending redward of the YSOs. Such red J − K_s colors, combined with relatively blue [4.5]–[24] colors, would be consistent with our hypothesis that many of these probable contaminants are AGB stars.

Zoom In Zoom Out Reset image size

Spectral index in the infrared, defined as

$$\begin{eqnarray}&&\alpha =\displaystyle \frac{d\mathrm{log}(\lambda {f}_{\lambda })}{d\mathrm{log}\lambda },\end{eqnarray} \tag{ 6 }$$

is frequently used to assess the evolutionary stages of YSOs (e.g., Lada 1987; Andre & Montmerle 1994; Evans et al. 2009a; Rebull et al. 2014). However, the value of α depends on what spectral range is used, with the largest available range typically being favored by most studies. Furthermore, the calculation of spectral index may also be affected by reddening. To estimate α values that are minimally affected by reddening, we use the wavelength range from 4.5 to 24 μm, since interstellar extinction in these bands is smaller than at shorter wavelengths and the reddening curve is flatter (Indebetouw et al. 2005; McClure 2009; Xue et al. 2016). For these bands:

$$\begin{eqnarray}&&{\alpha }_{[4.5]-[24]}\approx 0.55([4.5]-[24])-2.94\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{\alpha }_{[4.5]-W4}\approx 0.58([4.5]-W4)-2.92\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{\alpha }_{[4.5]-[8.0]}\approx 1.64([4.5]-[8])-2.82.\end{eqnarray} \tag{ 9 }$$

Where available, we prefer the α estimate based [4.5]–[24], followed by [4.5] − W4, and finally [4.5]–[8.0]. For YSOs suspected of having strong silicate absorption or PAH emission (Sections 5.6–5.7), we do not use the [4.5]–[8.0] color to estimate YSO class, because either feature could affect the 8.0 μm band.

Figure 15 shows the distribution of spectral indices calculated for candidate YSOs. Based on these estimates, there are 15,943 Class I (α > 0.3), 23,810 flat spectrum (0.3 ≤ α < −0.3), 59,949 Class II (−0.3 ≤ α < −1.6), and 5352 Class III (α ≤ 1.6) YSOs, using the α boundaries from Greene et al. (1994). In addition, there are 12,392 candidate YSOs with uncertain classes due to missing photometry. This classification scheme roughly reflects the YSO evolutionary sequence from deeply embedded sources with massive envelopes (Class I and flat spectrum) to stars with disks (Class II) and systems where the disk has mostly dispersed (Class III). However, viewing geometry may also affect the assigned YSO class: for example, a YSO that would otherwise be considered Class II may have a Class I SED if viewed at high inclination (Williams & Cieza 2011). Finally, we clarify that, even though some classification schemes regard Class III sources as having no IR excess (see Evans et al. 2009a), in our scheme Class III implies weak but detectable excess.

Zoom In Zoom Out Reset image size

Broad silicate dust absorption or emission features, centered at ∼9.7 and ∼18 μm, are frequently detected in the mid-IR spectra of YSOs (e.g., Furlan et al. 2006, 2008, 2011; Oliveira et al. 2010). The 9.7 μm feature overlaps the IRAC 8 μm band, so these features can affect YSO colors observed by IRAC.

In the color–color diagram shown in Figure 16 (left panel), a group of ∼2000 YSO candidates stand out due to their unusually blue [5.8]–[8.0] < 0 colors—these objects are flagged in Table 1. Given the lack of a red color in [5.8]–[8.0], the classification of these stars as YSO candidates was based mainly on their [3.6]–[4.5] ≳ 0.5 and [4.5]–[5.8] ≳ 0.5 colors, both of which tend to be redder than most of the other YSO candidates.

Zoom In Zoom Out Reset image size

Figure 17 shows three example SEDs that we have fit with YSO models from Robitaille (2017)—for each source, the 10 best-fitting convolved models are indicated by the gray lines. Robitaille (2017) include multiple configurations of disks and/or envelopes, so we used the simplest model forms capable of explaining the data: a star and disk model (sp-s-i ¹¹ ) for SPICY 75228; a star, disk, and envelope models with variable inner radius (spu-hmi) for SPICY 85135; and a star and disk model with variable inner radius (sp-h-i) for SPICY 99415. Although these fits are not all formally good, given the reported photometric uncertainties, they illustrate the range of SED morphologies that could produced the colors that we observe. Each case requires a strong silicate absorption feature at 9.7 μm to reproduce the lower 8.0 μm band emission. The best models also tend to also have nearly edge-on inclination to provide the high absorbing column density.

Zoom In Zoom Out Reset image size

Silicate absorption in YSO SEDs can come from the object itself or from foreground interstellar dust (van Breemen et al. 2011). We would expect a YSO with strong intrinsic silicate absorption to have a substantial disk or envelope that can produce the extinction, and Forbrich et al. (2010) find that YSOs with positive spectral indices are more likely to have strong silicate absorption. Figure 16 (right panel) shows our strong silicate absorption candidates on a plot of [3.6]–[4.5] versus [4.5]–[24]. The color [4.5]–[24] is a good indicator for SED spectral index that is not affected by silicate absorption. Most of the objects with possible silicate absorption have [4.5]–[24] > 5, higher than average for the YSOs, but ∼20 sources have [4.5]–[24] colors bluer than this. In the interstellar medium, the relation between the optical depth of the 9.7 μm feature and optical extinction is approximately τ_9.7 ∼ A_V /20 (Roche & Aitken 1984; Chiar et al. 2007; Shao et al. 2018). Although most stars in our sample would not have sufficiently high foreground extinction for the feature to become optically thick, this can be achieved along lines of sight that pass through dense molecular clouds or near the Galactic Center.

Another salient feature in the [4.5]–[5.8] versus [5.8]–[8.0] diagram is the finger-like structure at [5.8]–[8.0] ≈ 1.6. Most of the sources with these colors in the full IRAC catalogs were classified as probable contaminants, but a minority (∼490 objects) were classified as YSO candidates. These colors match those expected for sources dominated by PAH emission bands. For an astronomical PAH emission spectrum, the ratio of flux in the 5.8 μm band to the 8.0 μm band ranges from 0.31 to 0.41 (Draine & Li 2007, and references therein), corresponding to [5.8]–[8.0] = 1.6–1.9. There is little PAH emission in the 4.5 μm band, leading to a red [4.5]–[5.8] color. Candidates are flagged in Table 1 for strong PAH emission if they meet the criteria [4.5]–[5.8] > 2 and [5.8]–[8.0] > 1, which are based on the observed morphology of this feature in color space.

Although IR nebulosity in star-forming regions is dominated by PAH emission, inspection of the flagged YSO candidates suggests that most are valid point sources in all four IRAC bands, not spurious detection of nebular knots. For example, ∼90% of these sources have M = 2 detections in both the 5.8 and 8.0 μm bands, indicating reliable detections. We examined the images of a subset of these objects by eye and found that, even in cases with surrounding nebulosity, the sources themselves appeared to match the PSF. PAH emission may be intrinsic to massive YSOs with sufficiently high ultraviolet luminosities (e.g., Whitney et al. 2013). Spitzer IRS spectroscopy of massive YSOs has shown PAH emission to be nearly ubiquitous and correlated with YSO luminosity (Oliveira et al. 2013).

In the MYStIX IR-excess catalog that we used for training, Povich et al. (2013) aggressively filtered sources with PAH emission in order to avoid contamination by nebular PAH knots. To be included, they required sources to exhibit red K_s − [4.5] colors (avoiding bands with PAH emission), as would be expected for a massive YSO. This requirement will be reflected in our classifications via the random forest classifier. In the SPICY catalog, YSO candidates with possible PAH emission have median K_s − [4.5] = 3.3, as compared with a median for the entire sample of 2.0.

Figure 16 (right panel) shows the 24 μm emission for these objects. The YSO candidates with possible PAH emission have [4.5]–[24] colors ranging from 5 to 10, much higher than the average for YSOs. This result is consistent with these objects being massive YSOs.

Below, we demonstrate an example application for these cutouts, using a simple unsupervised image clustering strategy to characterize environments in which the YSOs candidates are found.

We avoid clustering in the pixel space because it is not invariant to image translations and rotations, which are properties that any proper content-based image clustering solution should have. Two candidate transforms that can introduce these properties via the power spectrum are wavelets (as used for a similar application in Krone-Martins et al. (2019)) and Fourier transforms (e.g., Kauppinen et al. 1995; van der Schaaf & van Hateren 1996); here, we adopt the latter. This is partially motivated because the Fourier power spectra is linked to the turbulent properties of the star formation medium (e.g., Elmegreen & Scalo 2004), revealing signatures of different physical phenomena.

We first compute the 2D Fourier power spectra of each cutout in each IRAC band. We then compute 1D radially medianized power spectra from each of the original 2D power spectra and concatenate these 1D power spectra to form a vector for each YSO candidate. Next, we organize the vectors of all environments into a single matrix and perform principal component analysis (PCA; Pearson 1901; Hotelling 1933), from which we select the most relevant dimensions (see also Ishida & de Souza (2013) and de Souza et al. (2014) for PCA variants), which acts as feature compression (see, e.g., Sasdelli et al. 2016). Finally, we model the distribution using a multivariate Gaussian mixture model (GMM; Pearson 1894; Scrucca et al. 2016; de Souza et al. 2017; Melchior & Goulding 2018) in the space defined by the first two principal components of the power spectra and the modes of the pixel values in each cutout, which we transform using an inverse sinh function. The distribution in this space is complex and requires many (25) Gaussian components, with model selection using the Bayesian information criterion (Schwarz 1978).

Visual inspection shows that the GMM components tend to correspond to three types of environments: those that are nebulosity-free (or have minimal nebulosity), mixed environments, and cloud-like environments. These are labeled Environments 1, 2, and 3 in Table 1. We also found outliers on the boundary of the distribution. Examination of the cutouts showed that the outliers correspond to severe image reconstruction errors and/or missing data in one or more bands and are located at the edges of the surveys.

A total of 66,539 stamps, or ∼57% of the valid stamps, were classified as cloud-like, while 32,790 stamps, or ∼28%, were classified as nebulosity-free. The mixed class and the outliers correspond to 15,462 and 2433 stamps, or ∼13% and ∼2%, respectively. These numbers indicate that most YSO candidates are indeed in cloud-dominated environments, as would be expected for YSOs in star-forming regions. However, the number of candidates in environments presenting diffuse nebulosity or even no detectable nebulosity is not negligible. There is a slight but statistically significant ¹⁴ correlation between spectral slope α and environment class, with stars in nebulosity-free environments having preferentially more negative spectral indices (indicating later evolutionary stages).

The cloud-like environments are most prevalent in the inner regions of the Galaxy, between approximately ℓ = 300° and 50° and ∣b∣ ≤ 1°. Cloud-like environments are also associated with large star-forming complexes outside this coordinate range, including some of the star-forming regions in Cygnus X.

The mixed environment is most prevalent further from the Galactic center (e.g., ℓ ≤ 310° or ℓ ≥ 30°). Some stellar associations include both cloud-like and mixed classes (e.g., the Carina Nebula).

The image cutouts with no or minimal nebulosity are found throughout the entirety of the survey, except within ∼1° of the Galactic Center. Stars in this environment are the most evenly distributed—but even among these stars, several clusters can be seen. For example, the Sco OB1 Association is made up of both cloud-like and nebulosity-free classes.

This simple application can certainly be significantly improved by adopting tailored methodologies and signal representations, for instance, curvelet transforms (e.g. Candès & Donoho 2000; Starck et al. 2003), to characterize the signal power contained in filamentary structures, and customized clustering methods. Moreover, a proper physical characterization of the YSO environment requires consideration of effects of, for example, the distances to the objects, accounting for the distinct physical scales probed, differences of PSFs in different IRAC bands and their impacts on the power spectra, projection effects of the nebulous matter in the plane of the sky, etc. However, as we show here, even a simple analysis has already revealed that, curiously, a significant fraction of the YSO candidates in the SPICY catalog do not seem to be lying in environments dominated by clouds.

To estimate the heliocentric distances (d_⊙) to each of the YSO groups, we employ a hierarchical Bayesian model (Hilbe et al. 2017) to account for the measurement errors in parallaxes as well as the presence of outliers. The use of robust statistics is particularly suitable given the non-negligible presence of unknown contaminants in each group. Normality assumptions are sensitive to noise and outliers, which may result in a biased estimate of the mean distance. Replacing a Gaussian likelihood with a t-distribution is a relatively easy fix. The t-distribution has an extra ν parameter called "degrees of freedom," which controls how closely the distribution resembles the normal distribution. Larger values ν > 30 essentially recover the normal distribution, while smaller values result in a distribution with heavier tails. This extra flexibility enables it to adapt to the extra noise in the data without introducing a bias in the underlying relationship.

The model formulation for the robust estimate is given below, where we define a t-likelihood for the observed ϖ and suitably vague priors on all the model parameters: uniform for d_⊙ over 25 kpc, and a gamma (Γ) prior (to ensure positivity) for ν.

$$\begin{eqnarray}\begin{array}{rcl}{\varpi }_{i} & \sim & {\mathscr{T}}(1/{d}_{\odot },{\sigma }_{{\varpi }_{i}}^{2},\nu ),\\ \nu & \sim & {\rm{\Gamma }}(2,0.1),\\ {d}_{\odot } & \sim & \mathrm{Uniform}(0,25),\\ i & = & 1...{n}_{\mathrm{Gaia}}.\end{array}\end{eqnarray} \tag{ 10 }$$

The index i runs over the members of each group n_Gaia with Gaia astrometric information. Although distance is constrained to be positive, our likelihood model permits the parallax measurements for individual stars, ϖ_i , to be either positive or negative. We evaluate the model using a Gibbs sampler, for which we use the jags ¹⁶ package (Plummer 2017) within the R language. We initiate three Markov Chains by starting the Gibbs sampler at different initial values. Initial burn-in phases were set to 5000 steps, followed by 5000 integration steps for each YSO group, which are sufficient to guarantee the convergence of each chain.

Table 2 provides group parallaxes and proper motions estimated from the posterior medians. Uncertainties are estimated from the mean absolute deviation (MAD) of the posterior (scaled to approximate 1σ uncertainties) and added in quadrature to the ±0.04 mas and ±0.07 mas yr⁻¹ spatially correlated systematic errors on DR2 zero points (Lindegren et al. 2018). Out of 406 groups, 402 have some Gaia astrometry, giving at least a rough estimate of parallax and proper motion. Of these, most groups include at least n_Gaia = 10 members having Gaia five-parameter astrometric solutions, enabling estimates that are more precise than those based on individual stars.

For each group, we show scatter plots of stellar proper motions, parallaxes, and positions (Figure 19), with the groups' mean parallaxes and proper motions indicated. In most cases, the stars form a single clump in ${\mu }_{{{\ell }}^{\star }}-{\mu }_{b}-\varpi$ space, suggesting that most of the group members are spatially and kinematically associated. In other cases (e.g., G77.8+1.0), multiple clumps are apparent, which may imply that distinct stellar groups with chance alignment have been merged by the HDBSCAN algorithm. In the example G77.8+1.0, the estimated properties correspond to the more distant but more numerous of the two groups. We visually inspected all groups in Figure 19 and have flagged those in Table 2 for which problems—such as the suggestion of multimodality, groups that appear dominated by field stars, or a single data point with too much leverage—could affect the interpretation of the Bayesian models. We also flag groups with n_Gaia < 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (406 panels)

Tarball of figure set images

The locations of these stellar groups in heliocentric Galactic XY coordinates are plotted in Figure 20. For this plot, we converted parallaxes to distance using an average −0.0523 mas Gaia DR2 zero-point correction estimated by Leung & Bovy (2019). Thus, $X={d}_{\odot }\,\cos (b)\cos ({\ell })$, and $Y={d}_{\odot }\,\cos (b)\sin ({\ell })$. The estimated centers of spiral arms from Reid et al. (2019) are also indicated. Fainter red points are groups with ϖ/σ_ϖ < 2 or groups that have been flagged. Few YSO groups are detected within 1 kpc, but the footprints of the GLIMPSE (and GLIMPSE extensions) surveys exclude many of the nearest star-forming regions, which are located more than several degrees above or below the Galactic midplane. The bulk of the YSO candidates for which we have accurate measurements have heliocentric distances that range from 1 to 3 kpc. There may be some bias in the distances to which we are sensitive, because this range resembles the range in the distances of objects in our training set. Nevertheless, there are groups that appear to lie beyond ∼3 kpc, but Gaia-based distances become more uncertain at this range.

Zoom In Zoom Out Reset image size

The YSO groups are not distributed smoothly within the Galaxy, but instead reveal Galactic structure. The survey areas intersect several spiral arms, and provide crucial information about Galactic structure at the boundary between Quadrants I and IV, where structures traced by v_lsr measurements of gas become degenerate. We discuss the relation of the stellar groups to the spiral arms below.

• Local (Orion) Arm: In Quadrants I and II, the Local Arm intersects both the SMOG field and the Cygnus X field. In SMOG, there is one group at the approximate distance of this arm. We find the YSO groups in Cygnus X to be spread linearly, spanning a factor of ∼2 in distance (∼1–2 kpc). This situation is consistent with looking down the length of the arm. Cygnus X is thought to lie at an end of a long molecular filament (Alves et al. 2020; Zucker et al. 2020) that contains multiple prominent star-forming regions (e.g., Orion, Taurus, the North America Nebula). The additional length of Cygnus X may increase the total length of this structure by 50%. A similar result is reported by Xu et al. (2016). In Quadrant IV, a chain of groups is located at X ∼ 0 and extends from ∼1 to ∼3 kpc in distance toward the constellation Vela. The orientation of this chain suggests that these groups could connect with the Local Arm.
• Sagittarius–Carina Arm: Numerous star-forming regions can be found within 20° of the direction of the Galactic center, including some famous regions like the Trifid Nebula, the Lagoon Nebula, and NGC 6334. Some of these groups (including the aforementioned famous regions) form a coherent, chevron-shaped structure that has a vertex pointing toward us at a distance of ∼1.0 kpc and edges that extend away with lengths of ∼1 kpc. This collection of star-forming regions is at the approximate distance of the Sagittarius–Carina Arm from Reid et al. (2019), implying that the arm is an active site of star formation activity. However, the angles of the edges of the chevron are inconsistent with the angle of the arm predicted by Reid et al. (2019). The linear structure making up the edges of the chevron cannot be a result of the "Fingers of God" effect, because it is not oriented along our line of sight. We find relatively little sign of YSO groups associated with the Sagittarius–Carina Arm at Galactic longitudes beyond ℓ > 30° in Quadrant I or between 300 < ℓ < 330° in Quadrant IV.
• Scutum–Centaurus Arm: This arm is less clearly delineated by stellar groups than the others, possibly owing to large distance uncertainties at the distance of this arm. However, there is an increase in the density of groups near this arm in Quadrants I and IV.
• Perseus Arm: This arm is intersected by the SMOG field, and three groups have distance estimates consistent with the center of this arm from Reid et al. (2019). The large distance of this arm may decrease our sensitivity to YSOs associated with it.
• Inter-arm: There are multiple stellar groups that appear to be located between the spiral arms from Reid et al. (2019). For example, within ∼10° of the Galactic Center, many groups appear located between the Sagittarius–Carina and Scutum–Centaurus arms. In Quadrant I, several groups are located between the Local Arm and the Sagittarius–Carina Arm.

The procedure for calculation of mean proper motions for the groups is similar to the calculation of heliocentric distances, using a weakly Gaussian prior for ${\mu }_{{\alpha }^{\star },0}$ and μ_δ,0 instead.

Figure 21 displays the proper motions in Galactic longitude and latitude as function of ℓ. The expected distribution for stars in circular Galactic orbits would be governed by Galactic rotation, parameterized by the Oort (1927) constants, and effects from solar motion, which are distance-dependent (Bovy 2017). On the ${\mu }_{{{\ell }}^{\star }}$ versus ℓ diagram, to first order in the Galactic plane, this would be a sinusoid with period 180° for Galactic rotation added to a sinusoid with period 360° for solar motion. This overall structure appears to dominate the ${\mu }_{{{\ell }}^{\star }}$ versus ℓ plot for our YSO groups, but there are hints of deviations that we will discuss in more detail in a subsequent paper.

Zoom In Zoom Out Reset image size

On the μ_b versus ℓ diagram, several structures can be seen. For example, between ℓ ∼ 5°–25°, there is a diagonal chain of groups in position–proper-motion space. These groups correspond to one of the edges of the chevron-like structure detected in the XY diagram. The dispersion in μ_b is slightly higher around ℓ ≈ 75° (Cygnus X) and around ℓ ≈ 260° (Vela)—both of which correspond to extended structures along the line of sight.

YSO candidates of all SED classes are spatially clustered, but the classes corresponding to earlier evolutionary stages tend to be more strongly clustered, as can be seen in the section of the Galactic midplane in Figure 22 (left panel). In this region near the young cluster NGC 6823, Class I sources mostly lie within the densest groups, but the other classes are comparatively more distributed. The K function (Ripley 1976) can be used to quantitatively compare the relative strength of clustering for these populations. We used the spatstat package (Baddeley 2017) to estimate K as a function of angular separation r for Class I, flat spectrum, Class II, and Class III objects with 95% confidence intervals estimate using the bootstrap method of Loh (2008). On this log–log plot (Figure 22, right panel), at angular separations of several arcminutes, the slope for Class I YSOs is significantly flatter than for flat spectrum YSOs, which is also significantly flatter than those for Class II and III YSOs, implying that the earlier stages are more clustered. This finding agrees with numerous other examinations of the spatial distribution of sources by YSO class (e.g., Sung et al. 2009; Samal et al. 2010; Buckner et al. 2020). We note that even some candidate Class I YSOs appear isolated: for example, ∼100 of these objects (<1% of the Class I YSOs) are separated from their nearest neighbors in our catalog by more than 10'.

Zoom In Zoom Out Reset image size

Figure 23 shows the smoothed distributions of sources of various classes in both Galactic longitude and latitude. In longitude, the normalized distributions of YSOs of all SED classes are similar, whereas in latitude, the distributions of earlier classes (e.g., Class I and flat spectrum) are more strongly concentrated near the midplane than those of the later classes (e.g., Class II and III). This may be a result of the dispersal of YSOs, if stars are born in regions nearest the midplane and then drift away. For example, a YSO traveling at a tangential velocity of ∼2 km s⁻¹ at a distance of ∼2 kpc could travel 025 from its point of origin in ∼5 Myr. This is enough to flatten the distribution of b shown in the figure but not the distribution of ℓ.

Zoom In Zoom Out Reset image size

The sources with strong silicate absorption are more concentrated toward the Galactic center than other YSOs. This could be an effect of the higher interstellar dust column densities in this direction. The YSOs with strong PAH emission also appear to be preferentially concentrated toward the inner Galaxy, but the peak of the distribution appears to be in star-forming regions around ℓ ∼ 330°.

Certain evolved stars, including dusty red giants, AGB stars, post-AGB stars, and red super giants (RSGs), have red IR SEDs due to dusty stellar winds (e.g., Marengo et al. 1997, 1999; Groenewegen 2012; Chun et al. 2015; Suh 2020), making such objects a significant category of contaminant in infrared YSO catalogs (e.g., Robitaille et al. 2008; Povich et al. 2013). Reiter et al. (2015) present a sample of AGB stars (including O-rich, S-rich, and C-rich stars) and RSG stars with JHK and IRAC photometry. In the near-IR, the J − K colors of this sample (J − K ≳ 0.9) are consistent with the group of probable contaminants on the J versus J − K_s diagram that are brighter and redder than the typical YSOs (Figure 9, left panel). In IRAC color space, the distribution of the AGB+RGS sources partially overlaps the distribution of YSO candidates. However, most of them have IRAC colors (e.g., [3.6] − [4.5] ≲ 0.5, [4.5] − [8.0] ≲ 1) bluer than the typical colors of YSOs but similar in color to the bright sources classified as probable contaminants (e.g., Figure 9, right panel).

Extragalactic sources, including active and star-forming galaxies, can have mid-IR colors that mimic the IR excesses of YSOs (e.g., Stern et al. 2005; Jarrett et al. 2011). These sources can contribute a significant number of possible contaminants in some YSO searches (e.g., Harvey et al. 2007; Gutermuth et al. 2008). However, in the GLIMPSE survey, extragalactic contamination is expected to be lower, owing to shallower IRAC observations and high extinction near the Galactic midplane (e.g., Kang et al. 2009). Jarrett et al. (2011) provide a sample of IRAC sources in fields dominated by extragalactic sources. The galaxies in these fields tend to have [3.6]–[4.5] ≈ 0–1.25 and [4.5]–[8.0] ≈ 1–4. This distribution more closely resembles the distribution of probable contaminants in our sample than it does our YSO candidates (Figure 11, lower right). Their extragalactic sources mostly have [3.6] ≳ 14 mag, which is fainter than most of our candidate YSOs.

For the random forest classifier, the sample of non-YSOs in our training set is equal in importance to the sample of YSOs. As with the labeled YSOs, the labeled field objects are obtained from the set of NIR+IRAC sources that could not be fit by a reddened stellar photosphere (Section 3.1). Thus, they also represent objects with red IR colors.

We generate our list of "non-YSOs" using both sources within the boundaries of the MYStIX star-forming regions that were classified as non-YSOs and field stars in regions of the Galaxy where there is no star formation activity. For the first category, we include all IRAC sources that lie within the MYStIX fields that show no indication of youth—requiring both rejection as a YSO based on its SED and nondetection of X-ray emission.

We also include stars from rectangular regions of the sky in areas where there is no evidence for star formation or the presence of YSOs. These regions, listed in Table 3, have been chosen to sample stars along different Galactic latitudes and longitudes and along lines of sight with different levels of extinction. Given the high numbers of stars within these fields, a random subsample of these stars are included in the training data for the classifier. Altogether, there are 14,019 labeled field objects for the 2MASS+IRAC sample, 5320 for UKIDSS+IRAC, and 4047 for VVV+IRAC.