The APOGEE-2 Survey of the Orion Star-forming Complex. II. Six-dimensional Structure

2.1.1. Data Products

2.1.2. Stellar Parameters

All the spectra are originally processed by the APOGEE Stellar Parameter and Chemical Abundances (ASPCAP; García Pérez et al. 2016) pipeline, providing various stellar properties, including the abundances (Holtzman et al. 2015). However, ASPCAP is primarily optimized for giants; given that the evolved YSOs targeted toward Orion are primarily dwarfs, this could introduce potential biases or systematic errors in the derived parameters. It is also a known issue that the uncertainties of various quantities (e.g., RVs) may be significantly underestimated, or that there may be a correlation in RVs for sources with T_eff < 3000 K (Nidever et al. 2015).

To account for these and other issues, we processed the spectra using the pipeline developed by Cottaar et al. (2014), which is better suited for the analysis of spectra of YSOs. It fits to each spectrum the effective temperature (T_eff), surface gravity ($\mathrm{log}\,g$), RV, rotational velocity ($v\,\sin \,i$), and veiling (r_H), using a synthetic grid and interpolating between the grid points. A Markov chain Monte Carlo (MCMC) simulation is performed at the end of the fit to measure the uncertainties in all the parameters and validated against the scatter in each parameter across multiple epochs of spectroscopic observations. Average parameters are determined for each star as a weighted average of the parameters measured for each epoch, as presented in Table 1.

Table 1.  Properties Derived from the APOGEE Spectra toward the Stars in Orion

α	δ	2MASS ID	Gaia DR2	Nepoch	Teff	$\mathrm{log}\,g$	$\overline{{RV}}$	${F}_{\mathrm{bol}}$	AV	Θ	Flaga	Group
(J2000)	(J2000)		ID		(K)	(dex)	(km s−1)	(10−10 erg s−1 cm−1)	(mag)	(μas)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)
85.1387100	−10.148064	2M05403328-1008530	3011876641301608832	1	4575 ± 46	3.360 ± 0.113	−12.379 ± 0.390	2.37 ± 0.10	2.40 ± 0.25	20.13 ± 0.59	0	field
85.5587769	−9.938928	2M05421410-0956201	3011907015310374528	1	4645 ± 31	4.326 ± 0.075	21.491 ± 0.528	6.51 ± 0.43	4.05 ± 0.20	32.39 ± 1.15	2	l1647-2
85.8625565	−9.993770	2M05432701-0959375	3011892137543646080	1	5388 ± 68	4.710 ± 0.092	20.863 ± 0.560	6.62 ± 1.49	2.05 ± 1.00	24.27 ± 2.81	2	l1647-2

Note.

^a0 = photosphere; 1 = photosphere with poor goodness of fit; 2 = photosphere up to H band.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

Table 2.  Astrometric Solutions from Gaia Utilized in the Analysis

α	δ	Gaia DR2	π	μα	μδ	Group
(J2000)	(J2000)		(mas)	(mas yr−1)	(mas yr−1)
(1)	(2)	(3)	(4)	(5)	(6)	(7)
85.2329562	−10.639077	3010342680846510208	2.262 ± 0.093	0.768 ± 0.170	−0.601 ± 0.158	l1647-3
87.9103488	−10.637480	3010909277228047744	2.385 ± 0.075	−0.358 ± 0.126	3.808 ± 0.111	field
86.0242576	−10.557282	3011793525094619008	2.109 ± 0.025	0.179 ± 0.045	−1.225 ± 0.043	l1647-2

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

To test the performance of the pipeline, we processed the APOGEE spectra of the stars observed toward the Pleiades with various synthetic grids, namely, the Coelho (Coelho et al. 2005), BT-Settl (Allard et al. 2012), and PHOENIX (Husser et al. 2013) spectral libraries, all of which were restricted just to the solar metallicities. A comparison was also made to the parameters previously obtained for the sources in the original IN-SYNC survey (Cottaar et al. 2014), which used the BT-Settl grid in spectral fitting. All the grids produced largely comparable solutions, although there were some systematic differences, most likely due to the underlying assumptions in the models (Figure 2). The most notable difference is in the absolute value of $\mathrm{log}\,g$: while the correlation between all models is largely linear, there is an offset of 0.5 dex between the solutions produced with PHOENIX and BT-Settl grids and an offset of 1.5 dex between the PHOENIX grid and the ASPCAP catalog. There are some differences in the absolute calibration of T_eff. The choice of grid has little to no effect on the determination of RV and $v\,\sin \,i$. We ultimately adopted the PHOENIX grid for our spectral analysis, as it produced the most self-consistent solutions and covers the widest range in T_eff (2300–15,000 K).

Zoom In Zoom Out Reset image size

While the pipeline does try to interpolate between grid points, it is less successful in some regimes than others (mainly with respect to T_eff), and particular values for the parameters are preferred to the exclusion of others. The most prominent zones of avoidance (nodding) for the PHOENIX grid occur at T_eff ∼ 3900, 5000, 5200, and 5700 K (Figure 3). A similar nodding effect is present in correlations with all grids, although it appears to be most frequent in the BT-Settl grid and weakest in the Coelho spectral library. Correlations with the PHOENIX grid also show a peculiar absence of sources at T_eff ∼ 3500 K, and to a lesser degree at T_eff ∼ 5000 K, although this gap does not follow the grid edges. The BT-Settl grid does not appear to show this gap (resulting in a dip deviating from the line of equity in Figure 2), and the presence or absence in the Coelho and ASPCAP grids is inconclusive owing to the limited temperature range. Comparison with optical surveys (e.g., Fang et al. 2018) suggests that most likely this 3500 K gap is due to the PHOENIX model producing NIR spectra that are not entirely representative of the real stars in this temperature regime.

Zoom In Zoom Out Reset image size

Veiling is a property that is only applicable to the stars that are accreting, and while a number of sources in the sample are disk-bearing accretors, a majority of them are not. However, when present, veiling will weaken spectral lines. Since the spectra were fit with only solar-metallicity models, the pipeline uses veiling to address any discrepancies in the absolute line depths. At low T_eff, the veiling parameter does systematically rise up to 0.3–0.5, even though most sources are expected to have veiling ∼0 (however, the accreting sources with high veiling can be identified above the systematic values). We correlated a small fraction of the stars (∼200 sources in Orion) with a larger PHOENIX grid that does include the effects of metallicity. The systematic effect in veiling at low T_eff for nonaccretors is significantly reduced to only ∼0.1. The correspondence in T_eff between models with and without metallicity variations is mostly one-to-one (although there is some scatter beyond 5000 K). Comparing the reduction where metallicity was and was not allowed to vary, the surface gravities are consistent for sources with $\mathrm{log}\,g$ > 4, but offsets grow for lower $\mathrm{log}\,g$. Below $\mathrm{log}\,g$ ∼ 3, the parameters derived assuming solar metallicity produced $\mathrm{log}\,g$ about 1 dex lower than those where metallicity is allowed to vary (primarily affecting the location of the red giant branch). The effects of nodding and the gaps described above remain unchanged.

The pipeline is capable of two modes of operation: fitting a single epoch of a given star at a time, or processing all the epochs to simultaneously fit a single T_eff, $\mathrm{log}\,g$, and $v\,\sin \,i$. However, with a large number of epochs, the latter mode produces increasingly unreliable solutions for the quantities that are allowed to vary between epochs, such as RV and veiling. We tested the pipeline on sources in the NGC 188 and Kepler 21 regions with 20+ epochs. Almost all sources had RVs of several epochs with signal-to-noise ratio (S/N) > 20 that were artificially scattered to values far outside the mean RV, compared to the solutions produced in the single-epoch fits. This effect was also present in sources with as few as three epochs. In extreme cases, an excessive noise from a single low-S/N spectrum may affect the invariant parameters as well. For this reason the parameters presented in this paper were derived by fitting each visit spectrum individually and then computing a weighted average of the visit-specific parameters, as the parameters from the individual visits have a much greater reliability. This is contrary to the approach taken in the IN-SYNC analyses, which adopted simultaneous fits that can result in spurious RVs. Because of this, identifications of spectroscopic binaries made on the basis of the original IN-SYNC multiepoch RVs may therefore be unreliable.

2.1.3. Uncertainties

The formal uncertainties σ in all the parameters are calculated by $\sigma =3{\sigma }_{\mathrm{MCMC}}\sqrt{{\chi }_{f}^{2}}$, where σ_MCMC is the variance in the MCMC output and ${\chi }_{f}^{2}$ is the reduced ${\chi }^{2}$ of the best fit. These calculations are similar to the approach by Cottaar et al. (2014) and provides a good agreement to the scatter in the data. If multiple observations of the same source are available, weighted averages (and the weighted uncertainties in those averages) are computed for all parameters from spectra with S/N > 20. It should be noted that these σ incorporate only the uncertainties in the fit, and not the systematic differences between various models.

For all sources with more than four epochs, we compute a degree of scatter ξ = Δ/σ, where Δ is defined as the absolute difference between the weighted average of a parameter and the value of that parameter at a given epoch (Figure 4). The distribution of ξ can be approximated by a Lorentzian 1/(γ(1 + ξ/γ)²) with a width γ = 0.4, and it is comparable for all parameters (the ξ_RV distribution does have a long tail, however, largely due to the presence of spectroscopic binaries). It is also largely invariant from the number of epochs available, or from which cluster a source originated. Approximately 80% of the scatter is within 1σ of the weighted average, ∼90% is within 2σ, and ∼95% is within 3σ. The median single-epoch uncertainties in all the parameters are ${\bar{\sigma }}_{\mathrm{RV}}$ = 0.4 km s⁻¹ , ${\bar{\sigma }}_{\mathrm{log}g}$ = 0.1 dex, ${\bar{\sigma }}_{{T}_{\mathrm{eff}}}$ = 50 K, ${\bar{\sigma }}_{v\sin i}$ = 1.7 km s⁻¹, and ${\bar{\sigma }}_{\mathrm{veil}}$ = 0.035. These uncertainties should be considered as a lower limit since they do not take into account potential added systematic offsets in the parameters.

Zoom In Zoom Out Reset image size

There does not appear to be a significant variation in the distribution of ξ of any parameters as a function of temperature or surface gravity at T_eff < 6250 K and $\mathrm{log}\,g$ > 1 (Figure 4, although RVs of giants with $\mathrm{log}\,g$ < 3.5 may be affected by jitter). Since stars with high T_eff have only a few spectral features that do not change significantly between different templates, σ_RV, σ_vsini, σ_veil, and σ_logg appear to be underestimated by a factor of 2, and ${\sigma }_{{T}_{\mathrm{eff}}}$ by a factor of 8 in sources with T_eff > 6250 K. We correct the uncertainties by the corresponding factors. The uncertainties for parameters with $\mathrm{log}\,g$ < 1 are somewhat more difficult to correct owing to a very low number of the affected sources; however, all of them are on the red giant branch, for which more precise fits are available in the ASPCAP catalog.

2.1.4. Stellar Radius via Spectral Energy Distributions (SEDs)

The observed APOGEE spectra have been fitted by the method discussed by Stassun & Torres (2016) and Stassun et al. (2017), in which a star's angular radius, Θ, can be determined empirically through the stellar bolometric flux, ${F}_{\mathrm{bol}}$, and effective temperature, T_eff, according to

$$\begin{eqnarray}&&{\rm{\Theta }}={({F}_{\mathrm{bol}}/{\sigma }_{\mathrm{SB}}{T}_{\mathrm{eff}}^{4})}^{1/2},\end{eqnarray} \tag{ 1 }$$

where σ_SB is the Stefan–Boltzmann constant.

${F}_{\mathrm{bol}}$ is determined empirically by fitting stellar atmosphere models to the star's observed SED, assembled from archival broadband photometry over as large a span of wavelength as possible, preferably from the ultraviolet to the mid-infrared. T_eff is measured from the APOGEE spectra, and thus the determination of ${F}_{\mathrm{bol}}$ from the SED involves only the extinction (A_V) and an overall normalization as the two free parameters.

We adopt the broadband photometry available in the literature spanning the wavelength range of approximately 0.15–22 μm. As demonstrated in Stassun et al. (2017), with this wavelength coverage for the constructed SEDs, the resulting ${F}_{\mathrm{bol}}$ are generally determined with an accuracy of a few percent. For stars with T_eff uncertainties of ≲1%, the ${F}_{\mathrm{bol}}$ uncertainty is dominated by the SED goodness of fit. With the exception of a few outliers, it was shown that one can achieve an uncertainty on ${F}_{\mathrm{bol}}$ of at most 6% for ${\chi }_{\nu }^{2}\leqslant 10$, with 95% of the sample having an ${F}_{\mathrm{bol}}$ uncertainty of less than 5%.

Here we apply this methodology to stars for which the estimated veiling in the H band was less than 0.5, suggesting little to no active accretion. Most of the stars should therefore be diskless, and the SED fit should be representative of the bare stellar photosphere. However, a number of the stars were found to exhibit significant excess emission in the WISE infrared bands, and in some cases in the Two Micron All Sky Survey (2MASS) K_S band as well. Therefore, we ran the SED fitting procedure a second time on these stars, now excluding the 2MASS K_S and the WISE photometry. Out of 8991 stars observed toward Orion, 6850 could be fit effectively as pure photospheres with the full wavelength range up to 22 μm, 298 had somewhat poorer goodness of fit, and 663 can be fit as photospheres up to H band, excluding longer-wavelength photometry. The remaining 1180 sources could not be fit owing to either high veiling or large infrared excess.

2.1.5. Sample

The sources observed toward the Orion Complex with APOGEE had substantial contamination from the foreground and background field stars. To minimize this contamination in the analysis throughout the paper, we excluded sources with RV < 10 km s⁻¹ or RV > 40 km s⁻¹, or those that fail the cuts in parallax, PMs, or color imposed in Section 2.2, if the data from Gaia are available (Figure 5).

Zoom In Zoom Out Reset image size

Gaia has recently produced its second data release (Gaia Collaboration et al. 2018), which contains astrometric solutions of parallax and PMs for 1.3 billion stars with G < 21 mag, significantly improving on the precision of previously available estimates of these parameters from Gaia DR1 and Hipparcos.

Previously, parallax and PMs have been measured for a number of young stars with the VLBA (Melis et al. 2014; Kounkel et al. 2017b; Ortiz-León et al. 2017a, 2017b; Dzib et al. 2018; Galli et al. 2018). The precision of these observations is comparable to that of Gaia DR2. Currently, a direct comparison of the measured distances between two surveys can be performed for 55 stars.

In general, there does appear to be good agreement between the measurements for single stars, as well as spectroscopic and long-period binaries (Figure 6), i.e., the sources that do not strongly deviate from an approximation of a linear PM. The astrometric solutions between the two surveys can be described with

$$\begin{eqnarray*}&&{\pi }_{{\text{}}{Gaia}}=(0.9947\pm 0.0066){\pi }_{\mathrm{VLBA}}-0.073\pm 0.034\end{eqnarray*}$$

$$\begin{eqnarray*}&&{\mu }_{{\alpha }_{{\text{}}{Gaia}}}=(0.9964\pm 0.0015){\mu }_{{\alpha }_{\mathrm{VLBA}}}-0.030\pm 0.073\end{eqnarray*}$$

$$\begin{eqnarray*}&&{\mu }_{{\delta }_{{\text{}}{Gaia}}}=(0.9960\pm 0.0020){\mu }_{{\delta }_{\mathrm{VLBA}}}-0.652\pm 0.034.\end{eqnarray*}$$

The systematic offset between parallax measurements is consistent with the zero-point offset of −0.03 mas reported by Lindegren et al. (2018). The large offset in μ_δ is driven by a few sources that may be affected by multiplicity that has not been apparent in the time frame of the VLBA-only observations; this offset is consistent with 0 for the remaining sources. As Gaia DR2 still does not have a prescription for astrometric binaries, the systems with a period of a few to a few dozen years present a considerable source of uncertainty that significantly increases the scatter in all parameters. In future data releases, this is expected to be improved.

Zoom In Zoom Out Reset image size

To analyze the structure within the Orion Complex, we selected the sources in Gaia DR2 with a position on the sky of 75° < α < 90° and −11° < δ < 15°, parallax 2 mas < π < 5 mas, and PMs −4 mas yr⁻¹ < μ_α < 4 mas yr⁻¹ and −4 mas yr⁻¹ < μ_α < 4 mas yr⁻¹ (Table 2). The ranges were chosen from the visual examination of the data to include all the overdensities that are associated with Orion. We excluded the astrometric measurements for the sources that had σ_π > 0.1 mas or σ_μ > 0.2 mas in order to retain only the sources with precise measurements. Furthermore, we discarded the sources that do not satisfy

$$\begin{eqnarray*}&&{M}_{G}\lt 2.46\times [{G}_{B}-{G}_{R}]+2.76;0.3\lt [{G}_{B}-{G}_{R}]\lt 1.8\end{eqnarray*}$$

$$\begin{eqnarray*}&&{M}_{G}\lt 2.8\times [{G}_{B}-{G}_{R}]+2.16;1.8\lt [{G}_{B}-{G}_{R}]\end{eqnarray*}$$

$$\begin{eqnarray*}&&{M}_{G}\gt 2.14\times [{G}_{B}-{G}_{R}]-0.57;0.5\lt [{G}_{B}-{G}_{R}]\lt 1.2\end{eqnarray*}$$

$$\begin{eqnarray*}&&{M}_{G}\gt 1.11\times [{G}_{B}-{G}_{R}]+0.66;1.2\lt [{G}_{B}-{G}_{R}]\lt 3,\end{eqnarray*}$$

where M_G = G + 5 − log(1000/π) in order to minimize the contamination from the main-sequence and red giant stars (Figure 7).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Throughout the paper, all the calculations are done in the parallax space. Conversion to distances occurs only when we report the averages and to obtain the stellar ages. When converting from parallax to the distance space, we correct for 0.03 mas offset and assume a systematic error of 0.02 mas (Lindegren et al. 2018).

We estimate the stellar ages using several different methods. One set of estimates is purely photometric. For YSOs without a reliable parallax measurement, we calculate M_G by assigning the average distance of their corresponding group (See Section 4). These calculations of M_G ignore the effect of extinction: outside of the youngest regions that are still associated with the molecular gas, representing most of the footprint of the Orion Complex, the extinction is expected to be relatively small (A_V ∼ 0). We then interpolate M_G and [G_B − G_R] over the PARSEC-COLIBRI grid of isochrones (Marigo et al. 2017). Despite the latest bandpass definitions, there appears to be a systematic offset between the data and the grid; to correct for it, we offset the grid's M_G by +0.2 mag and [G_B − G_R] by +0.15 mag. We note that the absolute calibration of the resulting ages may be imprecise, but the relative ages across the Complex should be largely consistent outside of the regions strongly affected by extinction. We refer to the ages derived with this technique as Age_CMD.

In addition to the above calculations, we also obtained ages for the sources observed with APOGEE, working with the T_eff obtained from the spectra (see Section 2.1.2) and ${L}_{\mathrm{bol}}$. We work with the extinction values computed as described in Section 2.1.4, as well as those estimated comparing the observed G − J and G_BP − G_RP colors from Gaia DR2 and 2MASS catalogs with the intrinsic colors obtained by interpolating T_eff in a modified version of Table A5 from Kenyon & Hartmann (1995) to allow for the new bandpasses and then transforming the color excesses to A_V using the coefficients C_G = 0.9145, ${C}_{{G}_{\mathrm{BP}}}=1.039$, and ${C}_{{G}_{\mathrm{RP}}}$ = 0.601 and a relation of ${A}_{V}=2.283{(({G}_{\mathrm{BP}}-{G}_{\mathrm{RP}})-({G}_{\mathrm{BP}}-{G}_{\mathrm{RP}})}_{o})$. With these extinctions and using the parallaxes from Gaia DR2, we estimated the absolute M_G and M_H magnitudes, which were converted to ${L}_{\mathrm{bol}}$ by applying the bolometric corrections from Kenyon & Hartmann (1995). To obtain the ages, we interpolated ${L}_{\mathrm{bol}}$ and T_eff into the PARSEC-COLIBRI isochrones, as explained in Suárez et al. (2017). When possible, ages obtained from both bands were averaged, and we refer to them as Age_HR. In contrast to Age_CMD, because of the dependence on T_eff from APOGEE, Age_HR values are not available for the entirety of stars that are identified as members of the Orion Complex. Where they are available, however, they are somewhat more reliable for the regions that are still associated with the molecular gas.

Finally, for the sources observed with APOGEE, we also interpolate the spectroscopically determined T_eff and $\mathrm{log}\,g$ values directly to estimate the ages. It is difficult to interpolate over these parameters with the PARSEC grid owing to strong nonlinear offsets between the grid and the data. But the relative distribution of ages for the various populations that we can infer from the spectroscopic parameters is consistent with the photometric distribution.

We do not report age measurements for the individual stars owing to the multitude of the utilized methods and the scatter in the measurements, which we plan to study and minimize in future work, but we discuss the population-averaged ages (both Age_CMD and Age_HR) in Section 4. In most cases, both are comparable, with the main difference originating from a somewhat different sample size.

Hierarchical clustering algorithms require one or more rules (or "constraints") to define where to "cut the tree" so that a small number of physically reasonable clusters are produced. These rules can be algorithmic, as with the choice of critical minimum cluster population in the widely used single-linkage clustering procedure of Gutermuth et al. (2009). But the constraints can be based on disciplinary knowledge external to the data set under study. We choose to set four constraints: maximum cluster size, maximum velocity range, an outlier rejection we call "reach," and the minimum number of stars. We set the values of these constraints based on a detailed simulation study described in Section 3.2. In statistical parlance, this is an example of "constrained" or "semi-supervised" clustering (Basu et al. 2008).

We identify structures using the 6D data (α, δ, π, v_α, v_δ, v_r), first standardizing them by scaling each dimension using the standard deviation of the values within it to avoid a mixture of units between positional and kinematical components and to give all dimensions an equal weight. Then, we compute the Euclidean pairwise distance matrix (L₂ norm) between all the sources using the IDL routine distance_measure. Groups that are formed with this metric are invariant to rotations and translations in the p-space (Duda et al. 2001).

Given that astrometric and spectroscopic catalogs do not have a complete agreement between the sources that are included, APOGEE-only sources have only three dimensions of data (α, δ, v_r), and Gaia-only sources have five dimensions (α, δ, π, v_α, v_δ). 6D distances involving these sources become undefined. To incorporate all the available data, we compute distance matrices on all three permutations (3D, 5D, and 6D) of the catalogs separately. All three matrices have the same size to simplify their merger later on even though certain elements of each matrix may not be defined in cases where data are incomplete. Nonetheless, even though all variables have been normalized, the measures of distance in each matrix are not compatible: those that were measured from the data with the larger number of the available dimensions are systematically larger (e.g., by $\sqrt{5/3}$). Additionally, the difference in source selections may further bias the measured distances. To correct for this, each matrix is then normalized by the 2nth-smallest element in the matrix, where n is the number of sources that have complete data in a given number of dimensions. This is an approximation of the median linking length of a constructed dendrogram for each permutation. The distance matrices are merged into a single distance matrix so that a unified clustering calculation can be performed, keeping the value from the permutation with the largest number of available dimensions for each pair of stars.

Then, a modified version of the IDL routine cluster_tree is used to identify structures in the data using the computed distance matrix, using the pairwise average linkage. Note that average linkage gives more compact and reliable clusters than the more commonly used single-linkage ("friends-of-friends") algorithm (Everitt et al. 2011). By default, it continues to run until every single source in a catalog is joined into a single group. To minimize excessive merging of unrelated structures or contamination from field stars, we imposed constraints that prevent a star from joining into a group if (1) the resulting group would be bigger than 4° in diameter (∼20–30 pc at the relevant distance range), (2) the absolute velocity difference (in any of the three dimensions) is bigger than 9 km s⁻¹, and (3) after a group has begun growing, it cannot accept any new stars if the linking length necessary to do so relative to the median linking length binding all stars to a cluster is greater than the threshold of 1.2 (cluster reach). At the end of the process, each star has a unique group assignment, and groups with fewer than 10 members are rejected. All these parameters are chosen empirically through the testing of performance on a synthetic data set, which is described below.

We applied the average linkage clustering algorithm on the synthetic population described in the Appendix. This population consists of a uniform distribution of field stars into which is embedded a randomly generated cluster with parameters such as distance, characteristic cluster velocity, age, and a number of stars that follow the normal distribution of a given velocity dispersion and the characteristic size. These parameters are varied in a manner that is representative of the ranges appropriate for the young stars toward Orion (Figure 8). While a single cluster cannot compare to the real population in terms of complexity, it can be used to test the performance of the algorithm in the various regimes, and the recovery properties of a single cluster are comparable to the recovery properties of multiple clusters in a single population.

For the purposes of clarity, the term "cluster" will refer to the generated population of YSOs, and the term "group" will refer to the output of the algorithm. To test the performance of the algorithm, we analyzed results for 1000 random clusters (each with a unique field contamination) processed by the clustering algorithm under a variety of input parameters (e.g., maximum group size, maximum group velocity). We evaluated the success of the clustering algorithm by computing the following for the resulting groups:

• 1.  
  False-positive fraction (FP) is defined as a number of field stars falsely identified as group members relative to the total number of grouped stars.
• 2.  
  True-positive fraction (TP) is defined as the number of true cluster members in a group relative to the total number of true cluster stars.
• 3.  
  If FP > 0.5, then the group is considered as fake. In a given set of permutations, a single TP and a single FP are calculated, adding the values for all the real groups (i.e., FP < 0.5) found.
• 4.  
  Split cluster fraction (SC) is defined as $\sum {N}_{\mathrm{real}}/{N}_{\det }-1$, where $\sum {N}_{\mathrm{real}}$ is the total number of real groups identified in all runs and N_det is the number of runs in which at least one group was detected (which may be <1000, as in some configurations of clusters no groups can be identified).
• 5.  
  Similarly, fake cluster fraction (FC) is defined as $\tfrac{\sum {N}_{\mathrm{fake}}}{{N}_{\det }}$, where $\sum {N}_{\mathrm{fake}}$ is the total number of fake groups identified in all runs.
• 6.  
  Missing cluster fraction (MC) refers to the fraction of runs without a single group detection relative to the total number of permutations, i.e., $1-\tfrac{{N}_{\det }}{1000}$.

Prior to determining the exact set of conditions that would be appropriate to use to truncate the branches of the cluster tree, we explored the parameter space, varying one property of the algorithm at a time, in order to consider the effect they had (Figure 9). The goal of this exercise was to minimize the contamination while maximizing the fraction of true cluster members recovered by the algorithm to determine the appropriate thresholds for all the parameters.

Zoom In Zoom Out Reset image size

Cluster reach (tolerance relative to the median linking length in a group to reject outliers) has relatively little effect on all of the parameters. As it relaxes, contamination increases as well (FP increases by 2% from varying the reach parameter from 1 to 2; FC increases by 4%), and clusters are more likely to be split into several groups (MC increase of 5% over the same range), while having small gains in TP (TP increases by 2.8%).

Varying a maximum allowed size of a recovered group has the largest effect on TP and MC: if it is too restricting, so as not to encompass the entire cluster, cluster members are more likely to be missed. On the other hand, contamination levels rise relatively modestly with more tolerant thresholds. Varying the maximum size from 2° to 4° increases FP by 2%, increases TP by 10%, and decreases SC by 8%. This parameter was fixed to 4°, which was the extent of the simulation, and which is almost large enough to encompass an entire molecular cloud in the Orion Complex.

The maximum velocity cut has a qualitatively similar effect to the maximum size; however, relaxing the threshold beyond 9 km s⁻¹ increases the contamination significantly, with FC rising by 14% at 12 km s⁻¹. Finally, the critical number of stars—defined as the minimum number of stars that can be considered as a group at the end of the output—also has a significant effect. Allowing smaller groups increases TP through producing smaller groups that split from the main cluster, while also exponentially increasing FC.

With the default parameters listed in Section 3.1 applied to the synthetic data set, TP = 79.5% ± 0.5% (27,648 out of 34,784 cluster members identified), FP = 6.20% ± 0.14% (1830 out of 29,524 grouped stars were from the field), SC = 27.8% ± 1.7% ($\sum {N}_{\mathrm{real}}=1247;$ N_det = 976), FC = 2.9% ± 0.5% ($\sum {N}_{\mathrm{fake}}$ = 28), and MC = 2.4% ± 0.5%.

Next, we tested the effect of the specific cluster properties on the recovery. Holding one cluster parameter at a time fixed to a specific value and allowing others to vary in the manner described in the Appendix, we generated 1000 clusters for each iteration. As the number of stars in a cluster becomes larger, the cluster is increasingly more likely to split into several smaller groups: clusters with 75 stars have produced multiple groups in 85% of iterations. While many of these splits occur along one dimension, some of the resultant groups are significant. While clusters were randomly generated, some correlations between the positional and kinematic components may occur, and multiple groups may portray this correlation. Sometimes, a split may occur just through a single dimension, and a closer examination of the output is necessary to confirm its significance. On the other hand, clusters with fewer than 30 stars are less likely to be recovered in the first place: almost 10% of clusters with 20 stars have not produced any groups.

Cluster distance does not have a significant effect on any of the tests, nor does the cluster age up to 15 Myr, as clusters older than that would require different photometric cuts, and neglecting any effects of that would increase the dispersion of the older clusters beyond the random generation. Compact clusters with small velocity dispersion are more likely to be recovered, with TP > 80% for the typical cluster radii less than 2.5 pc (with the Gaussian distribution of stars in a cluster, a total cluster size is >10 pc), and for velocity dispersion less than 2 km s⁻¹. Those with velocity dispersion larger than 1.5 km s⁻¹ are likely to run into the edges of the maximum allowed velocity range and be split into multiple groups, with the increase of SC from ∼20% up to 60% at 2.5 km s⁻¹. It should be noted, however, that in the Orion Complex the ONC is the most massive cluster, the velocity dispersion of which is ∼2.5 km s⁻¹ (Kounkel et al. 2016); in other regions dispersion on the order of 1 km s⁻¹ is not uncommon (e.g., Briceño et al. 2007; Nishimura et al. 2015).

Then, we considered the effect of the field contamination, by increasing the number of field stars by a multiple factor to what has been previously adopted in the Appendix. If the field population is underestimated in the generation of the synthetic population (or, if the photometric cuts are modified to include YSOs older than 15 Myr), then TP does decrease somewhat (dropping to ∼74% if the field population is 3 × 4.5 × 10⁴ stars), with a slight increase in FP (up to 11% at the same level), but the main difference is the significant rise in FC (almost every run produces one or more fake groups), which requires additional vetting.

Finally, a test was performed for the confusion in assignment of the stars in a multicluster population. Two clusters were placed at a random position in the sky in a field population covering 8 × 8 deg². Other properties were allowed to vary as above. A group was identified as being associated with a particular cluster if more than 2/3 of all stars identified in the group originated from that cluster. In a run consisting of 10,000 simulations, we recorded the fraction of stars from a neighboring cluster that made it into a wrong group, which was 0.6% ± 0.01%. We also tracked the fraction of the permutations in which at least one of the identified groups contained a mix of stars such that no clear cluster assignment could take place (with at least one-third of sources in a group originating from one cluster, and one-third from the second cluster). This fraction was somewhat higher at 2.04% ± 0.14%, where 60% of the affected groups were split off from the main groups associated with each cluster, 20% had confident grouping of one of the clusters and confusion regarding the assignment of the second one owing to the contamination from the first, and in 10% of the cases both clusters are merged into a single group. Unfortunately, due to a multitude of the dimensions in which parameters could vary, it is difficult to provide the specific conditions of the minimum distance between two stellar populations before confusion would occur.

In the observations of RV toward σ Ori, Jeffries et al. (2006) have observed two distinct velocity components separated by 7 km s⁻¹, one of which was interpreted to originate from the population of stars associated with σ Ori, whereas the other component is attributed to Ori OB1ab. We obtain similar results in the RV observations with APOGEE. The two components appear to extend northwest significantly beyond the area originally covered by Jeffries et al. (2006). The component centered at v_r ∼ 3–5 km s⁻¹ (Figure 14.1) traces the population of YSOs surrounding the belt stars, as well as the OB stars such as η Ori, 22 Ori, 25 Ori, and ψ² Ori. This population includes stars associated with Orion OB1a and OB1b. It stretches further south beyond the area covered by the APOGEE footprint toward Rigel, which includes a few of the Orion outlying clouds (Alcalá et al. 2008; Biazzo et al. 2014), and some authors recently referred to the region as Orion X (Bouy & Alves 2015; Zari et al. 2017). On the other hand, the component centered at v_r ∼ 13 km s⁻¹ (Figures 11, 14.2) consists of the σ Ori cluster, and it stretches diagonally toward an area east of 25 Ori, concluding northeast of ψ² Ori, largely excluding any other notably bright OB stars. The two populations are also found at different distances; this has also been recently found by Briceno et al. (2018). This demonstrates that σ Ori originated from a different cloud than the rest of the OB1ab regions. We extend the naming convention of the Orion A and B molecular clouds, and we will refer to this progenitor as the (former) Orion C molecular cloud. Collectively, we will also refer to the entire region encompassing OB1a, OB1b, and Ori X as Orion D. We note that the original definition of the OB1 subassociations is largely based on the area on the sky—not necessarily representative of the underlying structure. As two populations are largely superimposed onto each other, caution must be exercised in comparison with the literature, especially in the vicinity of OB1b. The motivation behind grouping it with Ori D is due to the tighter correlation of the foreground population with the belt stars, but many of the stars that have been associated with OB1b are also associated with the Orion C region.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (5 panels)

Tarball of figure set images

At the present day, the gas from Orion C has been almost completely dissipated. There does appear to be a clear age gradient in the stars along the filamentary structure from α = 85°, δ = −3° to α = 82°, δ = 3° (Figure 13). Group C-North is the oldest, having a wide distribution in ages around a median of ∼7.5 Myr; C-Central has an equally broad distribution, with a median of ∼5.5 Myr, but also a peak at ∼3 Myr; and the σ Ori cluster is the youngest population, with an age of 1.9 ± 1.6 Myr. Both Age_HR and Age_CMD show a very similar distribution. Correcting for systematics, we measure the weighted average distances of 406 ± 4 pc for σ Ori, 413 ± 4 pc for Ori C-C, and 416 ± 4 pc for Ori C-N. Orion C also appears to be moving uniformly radially; however, it is expanding along the plane of the sky, with σ Ori and Ori C-N moving away in opposite directions from Ori C-C.

On the other hand, it is difficult to conclusively determine whether all parts of the Orion D population share the same progenitor or not, as the population containing them is largely continuous and mostly dispersed with few spatial kinematical differences. (In the case of OB1a and OB1b, this contradicts what has been previously observed by Briceño et al. [2007]; however, the sample of stars that they used as representative members of OB1b actually originated from the Orion C population.) Orion D has a comparable size in the sky to the entirety of Orion A, B, and C put together, and in another few megayears it is possible that those A, B, and C regions may become largely indistinguishable from each other. On the other hand, it is possible that Orion D is all part of a single expanding population, and most of the PM vectors are consistent with the process of expansion, with them dispersing most likely as a result of the loss of the molecular gas (e.g., Tutukov 1978).

Orion D has very little molecular gas remaining. Projected in part onto Orion B, the tail end of D corresponds to the very diffuse gas that can be seen in Figure 15 that has a gradient in v_r from 6 to 3 km s⁻¹. L1622 also is most likely associated with Orion D (see Section 4.3). A few diffuse clouds are found in the southwest, namely, L1616, L1634, and IC 2118. While most of the population is largely dissipated, some concentrated clusters are still apparent, such as 25 Ori at 354 ± 3 pc (Age_CMD = 6.2 ± 2.3 Myr, Age_HR = 7.4 ± 2.0 Myr), consistent with Briceno et al. 2018), ψ² Ori at 347 ± 3 pc (Age_CMD = 5.5 ± 1.7 Myr, Age_HR = 4.9 ± 1.5 Myr), η Ori at 347 ± 3 pc (Age_CMD = 5.7 ± 1.4 Myr, Age_HR = 6.4 ± 1.9 Myr), OB1b at 357 ± 3 pc (Age_CMD = 3.7 ± 1.8 Myr, Age_HR = 4.1 ± 1.6 Myr), and L1616 at 360 ± 3 pc (Age_CMD = 5.0 ± 1.7 Myr). The southern portion does deviate in distance from the rest of Ori D—the population that is associated with IC 2118 is located at 291 ± 2 pc (Age_CMD = 7.9 ± 2.5 Myr).

Zoom In Zoom Out Reset image size

A small group of stars are found north of κ Ori with a distance similar to IC 2118 at 302 ± 2.5 pc. This group has not been cataloged previously, although a few stars we identify have been previously confirmed as YSOs by Alcalá et al. (2000). We will refer to this group as Orion Y. Further tests will be necessary to confirm the properties and membership of this group. This group is spatially discontinuous from Ori D and older (Age_CMD = 10.5 ± 2.5 Myr). However, it does have similar PMs; therefore, we group it together with Orion D on all the plots.

Pillitteri et al. (2017) have also observed a nearby group of stars surrounding V1818 Ori for which they estimated a distance of 270 pc. It is not clear whether it is associated with Orion Y, as we do not recover it, and nearly all sources identified as YSOs by Pillitteri et al. (2017) have π < 2 mas.

One group identified by the algorithm is kinematically peculiar. 25 Ori-2 (it was recently identified by Briceno et al. 2018, under the name of HR 1833) is kinematically distinct by almost 10 km s⁻¹ in v_r and 1.8 mas yr⁻¹ in v_α from the main group. It also appears to be significantly older than 25 Ori, with Age_CMD = 15.1 ± 3.4 Myr (Age_HR = 12.9 ± 2.8 Myr). Its kinematics are similar to those of Orion C; however, its distance is more consistent with the Orion D population. Since it is closer than Orion C but its v_r is redshifted, it is unlikely that these stars have originated in Orion C; instead, they may be part of an earlier wave of star formation in Orion D.

The clustering algorithm used in this work has a relatively poor performance when it is applied on very numerous and extended stellar populations with a large velocity dispersion, creating many largely artificial distinctions between separate regions of the same cluster. As such, the ONC creates a particular challenge to disentangle in terms of its true internal structure (Figure 14.3). The high degree of extinction toward the cluster exacerbates the problem further. The ONC has been the subject of a number of studies to analyze its structure using methods that may be more suited to such an environment, including the work done by Hillenbrand & Hartmann (1998), Kuhn et al. (2014), Megeath et al. (2016), and Hacar et al. (2016), using 2D and 3D data. Here we present the identified groups mainly as a method to trace the extremes of the stellar population in a 6D space, with the caveat that the groups themselves may not be physical in nature, particularly compared to the results observed in the other regions in this paper.

In previous RV studies toward the ONC, a population of stars were identified that are blueshifted to the molecular gas from which these stars have formed, although no obvious correlations have been found with any of the stellar properties or the position of those stars on the sky (Tobin et al. 2009; Da Rio et al. 2016; Kounkel et al. 2016). Unfortunately, stellar PM information cannot be compared directly to the kinematics of the molecular gas. With the addition of the astrometric solutions from Gaia, the mystery of the blueshifted population remains, showing no correlation with the distribution of distances or PMs when applied to the optical observations from Kounkel et al. (2016). However, it is notable that the APOGEE sample does show significantly better agreement between the stellar RVs and those of the molecular gas.

In general, regarding the 3D structure of Orion A, there is a good agreement with the model from Kounkel et al. (2017b). Using Gaia astrometry, the average parallax toward the ONC is 2.540 ± 0.001 mas; corrected for systematics, this results in a distance of 389 ± 3 pc. The southern end of L1641 is found at π = 2.364 ± 0.004 mas, or 417 ± 4 pc, with relatively smooth continuity between them, and L1647 is found at π = 2.227 ± 0.006 mas, or 443 ± 5 pc.

Similarly good agreement is found for the individual measurements of kinematics (however, it should be noted that there was an error in conversion from μ_α,δ to v_α,δ in Kounkel et al. 2017b), but now it is possible to reconstruct the full map of the PMs along the cloud. The motions in the ONC appear to be mostly random, with slight preference for expansion near the outer edges. PMs in L1641 are preferentially oriented perpendicular to the filament in the plane of the sky (most of the grouped sources are located near the northern edge of the molecular cloud, and their PMs are oriented toward the gas); in L1647, stars are primarily moving away from the main filament in the plane of the sky.

We estimate Age_HR = 1.6 ± 1.5 Myr (Age_CMD = 3.2 ± 2.2 Myr, which is an overestimate due to extinction). This is consistent with the estimates by Getman et al. (2014). Similar distributions are present in the northern and southern parts of L1641, with ages of 2.0 ± 1.5 Myr and 1.9 ± 1.4 Myr (Age_CMD = 2.4 ± 1.8 Myr and 2.1 ± 1.7 Myr), respectively, consistent with the measurements from Hsu et al. (2013). L1647 is the youngest region, with Age_CMD = 1.3 ± 1.3 and Age_HR = 1.9 ± 0.9 Myr.

There is one group toward the ONC that is kinematically different from the rest (ONC-22), as it is offset from the main population in v_α by 2.6 mas yr⁻¹. Unfortunately, no RV measurements are available, but it does not significantly deviate in other dimensions, though it may be found closer to the back of the cluster.

Particular interest has been given in the past to NGC 1980, a population that exhibits relatively little extinction compared to the rest of Orion A. Several studies have suggested that this region is more evolved compared to the ONC and that it may be located somewhat in the foreground (Alves & Bouy 2012; Pillitteri et al. 2013; Bouy et al. 2014). However, there has also been some doubt to these claims (Da Rio et al. 2016; Fang et al. 2017; Kounkel et al. 2017a). We do find a diffuse distribution of YSOs in the foreground as part of the Ori D structure, and some of it does coincide spatially with Orion A. It does not appear to be particularly concentrated near the NGC 1980 region, and only ∼8% of the total that is found toward it (∼30 stars) is likely to be associated with the foreground population, but it is possible that they may have contributed to the confusion surrounding the issue.

Confusion with Orion D toward Orion A is most notable toward the north of the cluster, near NGC 1981, which is another population that has been previously observed to be older than the rest of the ONC (e.g., Maia et al. 2010). There are a number of groups that are observed to have a distance more consistent with that of Orion D, at 357 ± 3 pc, with Age_CMD = 5.0 ± 1.8 Myr and Age_HR = 2.8 ± 1.7 Myr. Most of these groups do not have reliable v_r measurements, but those that do have a good agreement with that of the ONC—somewhat surprising considering the peculiar RV structure of Orion A north of δ < −5°, which is not representative of Orion D. While we associate these groups with Orion D, this assignment may not necessarily be accurate.

Previously Kounkel et al. (2017c) studied the RV structure of the Orion B molecular cloud. They identified a peculiar RV structure toward the NGC 2024 region that is largely asymmetric from the RV structure of the molecular gas. Although due to a high degree of extinction toward it, their optical sample contained relatively few sources. With the expanded sample of the APOGEE observations, the agreement between the kinematics of the stars and the molecular gas is significantly improved. The blueshifted clump can be resolved into a kinematically distinct population that is centered at v_r ∼ 6 km s⁻¹ (Figure 15), which appears to be associated with Ori D. On the other hand, the redshifted clump is traced by the sources that originate from σ Ori (see Section 4.1).

Toward the NGC 2068 cluster the agreement between the stellar and molecular kinematics remains good. There is some hint of a slight excess of stars somewhat redshifted relative to the gas, but it is difficult to conclusively state how significant it is.

We do recover a population that is located east of L1617, in the gas-free region close to the edge of the molecular cloud. We will refer to this population as the Orion B-North group. All three clusters are found at comparable distances of 404 ± 5 pc, 417 ± 5 pc, and 403 ± 4 pc for the Ori B-N group, NGC 2068, and NGC 2024, respectively (Figure 14.4). There is a deviation between the distances measured by Gaia and those measured by Kounkel et al. (2017b) toward these regions; this could be attributed to the small sample size and multiplicity. In terms of PMs, both Ori B-N and NGC 2068 appear to be moving toward NGC 2024.

L1622 is found outside of the footprint covered by APOGEE, and we do not recover it owing to a limited number of sources in Gaia associated with it because of the high extinction. From seven Class II stars that have been identified by Megeath et al. (2012) as YSOs in L1622, we measure average astrometric parameters of π = 2.871 ± 0.026 mas (d = 345 ± 5.5 pc), v_α = 5.343 ± 0.045 mas yr⁻¹, and v_α = 4.271 ± 0.040 mas yr⁻¹. Combined with the typical v_r = 1.17 km s⁻¹ (Kun et al. 2008), we conclude that L1622 is not associated with the Orion B molecular cloud, but may possibly have been related to Orion D. On the other hand, while L1622 and L1617 appear to be projected in a similar area of the sky, it is probable that the two are unrelated. While no parallaxes are available for the sources associated with the molecular gas in L1617, it is more likely to be a part of Orion B given its v_r (Reipurth et al. 2008).

We estimate Age_HR of 1.1 ± 1.0 Myr and 1.0 ± 0.5 Myr (Age_CMD of 1.9 ± 2.0 Myr and 1.0 ± 1.4 Myr) for NGC 2024 and NGC 2068, respectively. Ori B-N group has Age_CMD = 2.2 ± 1.0 Myr.

Most of the molecular gas toward λ Ori has already been dissipated by a supernova; however, a couple of clouds within the ionized bubble remain, namely, B30 in the northwest, and B35 in the southeast. YSOs observed toward this region are primarily located alongside a filamentary structure stretching between these clouds. Stellar RVs are largely consistent with the RVs of the molecular gas, although a direct comparison is no longer possible (Figure 16).

Zoom In Zoom Out Reset image size

In the cluster centered on λ Ori, though, there do appear to be two distinct kinematical components (Figure 14.5): the main one at v_r ∼ 12 km s⁻¹, and a secondary subcluster (λ Ori-2) at v_r ∼ 14 km s⁻¹. The separation is also apparent in PMs, but not in parallax (average distance to the cluster of 404 ± 4 pc) or in the spatial position on the sky. It appears that λ Ori-2 favors somewhat older ages than λ Ori-1, with Age_HR = 3.7 ± 1.0 Myr versus 5.1 ± 1.1 Myr (Age_CMD = 3.5 ± 0.9 Myr vs. 4.4 ± 1.4 Myr). A few groups (λ Ori-3 and λ Ori−4) also have somewhat peculiar PMs, although they lack a reliable number of RV measurements. YSOs near B30 (located at a distance of 396 ± 4 pc) are younger, 2.4 ± 1.3 Myr, with a near-uniform distribution of ages from 2 to 5 Myr, with the northernmost sources (δ > 12°) near the cloud having Age_CMD = 1.2 ± 0.8 Myr (Age_HR = 2.1 ± 0.8 Myr). These YSOs are clumped preferentially near the inner edge of the molecular cloud, with a common v_r ∼ 10 km s⁻¹. A similar distribution is found toward B35, with a typical age of 2.6 ± 1.3 Myr (distance of 397 ± 4 pc), with the easternmost sources becoming somewhat preferentially younger with an age of Age_CMD = 1.8 ± 1.5 Myr (Age_HR = 2.0 ± 1.0 Myr).

A supernova was produced near the cluster center a few megayears ago. Nearly all PMs larger than 1 mas yr⁻¹ relative to the average cluster motion are moving radially away from λ Ori, and the farther the sources are away from the center of the cluster, the faster they tend to move (Figure 17). This is consistent with a single trigger expansion, which has a travel time from the cluster center of ∼4.8 Myr. On the other hand, sources within 15 from the cluster appear to be mostly relaxed.

Zoom In Zoom Out Reset image size

Based on their age, B30 and B35 would have formed halfway from λ Ori to their current position and would not have been accelerated owing to the sudden mass loss of the molecular gas from the cluster. The fastest-moving stars have PMs of up to 6 km s⁻¹ relative to the average motion of the cluster. Considering that

• 1.  
  observationally, kinematics of YSOs usually closely mirror kinematics of the gas from which they form,
• 2.  
  it is easier for an expanding shock wave to influence the velocity of the molecular gas compared to the already-forming YSOs, and
• 3.  
  outlying stars are significantly younger than those in the cluster center,

it is possible that the formation of a number of those stars has been triggered by the supernova, similarly to the scenario described by Mathieu (2008). However, it should be noted that from simulations it appears to be difficult to distinguish between stars whose formation has been truly triggered and those that are in the process of formation regardless of any feedback (Dale et al. 2015). Further modeling of the cluster dynamics will be needed to demonstrate the role that triggering has played in the formation of the outlying stars.

There may be additional structure near the outskirts of λ Ori, which was suggested by Zari et al. (2017; e.g., L1588), although we do not recover it.