IN-SYNC. V. Stellar Kinematics and Dynamics in the Orion A Molecular Cloud

The observational study of stellar kinematics in young clusters provides critical clues on the mechanisms governing cluster formation and early evolution. The velocities of young stars, in comparison with gas kinematics, can reveal if and for how long the newly formed stellar population follows the initial/early gas flows, which may be due to turbulence, gravitational infall, cloud collisions, or some other triggering mechanism. The star cluster may be forming in a quasi-monolithic fashion, or the process may involve merging of different subclusters. These mechanisms may be revealed by a study of the kinematics, including kinematic substructure, of the young stars. Subsequent dynamical evolution may involve mixing of orbits, e.g., as the star cluster virializes as a gravitationally bound system. Alternatively, depending on the overall efficiency of formation from the natal gas clump and on the degree of gas expulsion, the stellar population may be or become unbound and thus expand, with fewer opportunities for stellar interactions. Depending on the rate of star formation, i.e., the efficiency per local free-fall time, these potential evolutionary pathways may be followed even while the cluster/clump system is still gas dominated. Another way to test different star formation scenarios is thus to measure the velocity dispersion of the cluster and compare it with the value expected from virial equilibrium given the gravitational potential due to the total mass (stars + gas) of the system.

Great effort has been made over the years to study the dynamical evolution of young clusters theoretically, through N-body simulations (e.g., Scally & Clarke 2002; Goodwin & Bastian 2006; Baumgardt & Kroupa 2007; Allison et al. 2009, 2010; Fellhauer et al. 2009; Parker & Meyer 2012; Parker et al. 2014; Farias et al. 2015, 2017). Results reveal a fairly complex picture, in that the dynamical state and evolution of the stellar population, the bound fraction, and the morphology depend on several factors: the star formation efficiency, the time of gas removal and its duration, and the initial configuration of hierarchical/substructured protoclusters.

Observationally, some works have analyzed the morphology of young clusters at different ages to constrain their dynamical evolution and initial conditions (e.g., Gutermuth et al. 2005; Schmeja & Klessen 2006; Da Rio et al. 2014b; Banerjee & Kroupa 2015; Jaehnig et al. 2015). However, observations of radial velocities in young clusters are a more powerful tool to probe the current dynamical state, search for and parameterize binary populations, and investigate spatially coherent velocity gradients or substructure that might give clues to the initial conditions or reveal multiple populations (Jeffries et al. 2006, 2014; Fűrész et al. 2008; Tobin et al. 2009; Cottaar et al. 2012; Rigliaco et al. 2016). In many cases, these radial velocities, v_r, have been obtained from optical high-resolution spectroscopy surveys that reach radial velocity precisions of ∼1 km s⁻¹ or better for individual sources. Such precision is needed to resolve the kinematics of nearby relatively low-mass regions of star formation (e.g., Orion), which have velocity dispersions of a few km s⁻¹. In young embedded clusters, optical spectroscopy can be limited by dust extinction and thus be unable to reach stellar members affected by high extinction.

In this context, with the INfrared Spectra of Young Nebulous Clusters (IN-SYNC) project, see Cottaar et al. (2014), which is an ancillary program of the SDSS-III, we obtained multiobject near-infrared (NIR) high-resolution spectra in three young clusters. This program, which used the Apache Point Observatory Galactic Evolution Experiment (APOGEE) spectrograph, allowed us to derive stellar parameters and to reach radial velocity precisions to ∼0.3 km s⁻¹ down to H = 12.5 mag. The first part of this survey covered the Perseus cloud, through their young clusters IC348 (Cottaar et al. 2014, 2015, hereafter Paper I and Paper III) and NGC 1333 (Foster et al. 2015, hereafter Paper II). We have found evidence for a supervirial stellar population in IC 348, whereas in NGC 1333 the stars are in close agreement with virial velocities, but the diffuse gas and dense gas have a significantly different velocity dispersion. Moreover, we found that stars in NGC 1333 show a similar velocity dispersion to the diffuse gas, whereas the dense cores appear in subvirial motions, unless their dynamics are significantly regulated by large-scale magnetic fields.

In a second part of IN-SYNC, we covered the Orion A molecular cloud. At the distance of ∼400 pc (Menten et al. 2007), this filamentary structure includes the Orion Nebula Cluster (ONC), the closest site of active massive star formation. Sparser populations are present to the north, in the upper sword, and NGC 1977, and to the south, with NGC 1980 and L1641. This population has been thoroughly studied at many wavelengths through photometry and spectroscopy (Hillenbrand 1997; Hillenbrand & Hartmann 1998; Fang et al. 2009, 2013; Da Rio et al. 2010, 2012; Robberto et al. 2010, 2013; Hsu et al. 2012, 2013). Memberships have been estimated through optical spectroscopy, the presence of IR excess emission (Megeath et al. 2012), and X-ray emission (Getman et al. 2005; Pillitteri et al. 2013). The population has a typical model-dependent mean age of $\sim 2\mbox{--}3\,\mathrm{Myr}$ throughout the region and in spite of the large extent of the structure (∼40 pc). The only exception is the population around NGC1980, which is somewhat older (Alves & Bouy 2012; Bouy et al. 2014; Da Rio et al. 2015). In the ONC, a large age spread compared to the local free-fall time is present, see Da Rio et al. (2010, 2014a, 2015), indicating a slow star formation process relative to the dynamical timescale of these star-forming clouds (Tan et al. 2006; Krumholz & Tan 2007).

The kinematic properties of the molecular gas in Orion A have been studied in great detail (e.g., Bally et al. 1987; Dame et al. 2001; Nishimura et al. 2015; Ishii et al. 2016) from tracers such as ¹²CO, ¹³CO, and C¹⁸O. The filament shows a large ∼10 km s⁻¹ change in radial velocity over its length, with the south end (the "tail") blueshifted compared to the north end. In addition, based on 3D dust mapping, Schlafly et al. (2015) found evidence that the tail is several tens of pc more distant than the northern end. Unless the Orion A cloud is stretching out along its orientation in plane of the sky—which is unknown because of the lack of accurate proper motion measurements—this would suggest that the Orion A filament is compressing along its length. As for stellar kinematics, Fűrész et al. (2008) and Tobin et al. (2009) conducted an optical survey measuring v_r in the northernmost third ($\sim 2^\circ$ long) of the region, centered on the ONC. In this region, they found that the velocities of the young stars follow the gas (¹³CO) velocity, although their distribution is slightly asymmetric with a broad blueshifted tail that is not present in the gas component. They also noted that the large-scale v_r gradient steepens north of the ONC, which had been interpreted as evidence for large-scale infall.

In our IN-SYNC survey, we obtained 4828 spectra of 2691 individual sources throughout the $\sim 6^\circ$ long region. The spatial distribution of these sources is shown in Figure 1, left panel, together with a ¹³CO map from Nishimura et al. (2015). In the previous paper of this series (Da Rio et al. 2015, hereafter Paper IV), we used these data to study the properties of the stellar population. We focused on the fitted stellar parameters (${T}_{\mathrm{eff}}$, A_V, $\mathrm{log}g$, v_r), positioned the sources in the Hertzsprung–Russel diagram (HRD), and assigned stellar ages.

Zoom In Zoom Out Reset image size

Our stellar sample was limited to $H\lt 12.5$ mag and assembled by prioritizing known members from a multitude of literature estimates; additional sources with unknown membership were then added if they were in the luminosity range reached by our survey. Our sample reaches masses as low as $M\sim 0.15$ ${M}_{\odot }$ for low values of A_V, which in turn covers values up to ∼20 mag with a mean A_V = 2 mag. The coverage of known members within the luminosity range of the survey is nearly complete throughout the region, with the exception of the central part of the ONC, where fiber collision constraints, combined with crowding, limited us to about 50% completeness. Stellar parameters are very accurate for cool stars (${T}_{\mathrm{eff}}\lt 5000$ K), with mean uncertainties in ${T}_{\mathrm{eff}}$ on the order of 50 K and mean and median errors in v_r of ∼1 and 0.3 km s⁻¹. Our precision worsens at early spectral types; however, because of the IMF, only a minority of members are in this range.

We found a clear confirmation of a genuine age spread throughout the region, from the correlation between HRD ages and surface gravity-derived ages, as well as an anticorrelation between ages and IR excess, and ages and extinction. We also confirmed that the non-embedded population of young stars around NGC1980 is older than the rest of the system, but suggested that it is not part of a well-separated cluster. This is because the kinematic properties of these sources are indistinguishable from those of younger more embedded stars at the same position. Based on stellar parameters as well as radial velocities, we found ∼400 new candidate members. Most of these turn out to be young diskless sources, in portions of the Orion A cloud that are not fully covered by previous studies to assign memberships.

In this paper we continue the analysis of our IN-SYNC data set, focusing on the kinematics and dynamics of the population.

Figure 1 shows a map of the Orion region as well as a position (decl.)–velocity diagram of our IN-SYNC targets in comparison with that of the molecular gas tracers. Figure 2 shows an animation of the full 3D (R.A., decl., and v_r) cube of the data. It is qualitatively evident that stars follow the global motions of the molecular gas throughout the region, as seen before limited to the north portion of the cloud ($-6^\circ \lt \delta \lt -4^\circ$, Tobin et al. 2009). Moreover, unsurprisingly, non-members tend to have velocities that are different from the bulk of the young population. Another feature is the fact that there appear to be more members with v_r lower than the peak value than stars with higher velocities, in line with the results from Tobin et al. (2009) that the stellar velocity distribution is asymmetric with a broader blueshifted wing.

To investigate this possibility, we construct positional maps of mean stellar velocities. First we exclude from our analysis all non-members; we also exclude sources with velocities strongly different from that of the bulk of the population (${v}_{r}\lt -20$ km s⁻¹ or ${v}_{r}\gt 35$ km s⁻¹) to exclude outliers that either have erroneous membership or velocity assignment, are binary systems, or ejected members. Then, for every angular position in $(\alpha ,\delta )$, we consider the sources within a circular aperture with radius of $20^{\prime}$, provided that there are at least 30 sources in the aperture, and we compute the weighted average of their v_r values, as well as the median and mode value. Last, we smooth each map with a kernel of $5^{\prime}$ of radius.

Figure 3 shows the derived maps, as well as the difference between stars and gas, color-coded according to v_r and density or flux. Visually, it appears that the difference stars-gas, which we referred to as ${\rm{\Delta }}{v}_{r}$, is on average more negative then positive, as the overall color leans toward the blue. In fact, the area-weighted mean value of ${\rm{\Delta }}{v}_{r}$ using the stellar mean velocity is −0.87 km s⁻¹, which decreases in absolute value to −0.44 km s⁻¹ , weighting the result on the stellar density. These values decrease to −0.66 and −0.25 km s⁻¹, respectively, when adopting the mode of the stellar velocities. Such difference is also evident from Figure 4, where we show the velocity distributions of stars and gas in different declination bins. Therefore, it appears that stars are systematically slightly blueshifted compared to the gas, and this is in part due to an asymmetry in the v_r distribution, as manifested by the difference between mode and mean v_r.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It is possible that this offset is not real and is simply due to zero-point offsets in v_r. However, by comparing our radial velocities with those from Tobin et al. (2009) for 550 matching sources, we find a difference in the values of 0.11 ± 0.05 km s⁻¹, much smaller than the offsets between stars and gas. Unless we both have a similar systematic error in v_r even though we used different instruments and techniques, this suggests that our IN-SYNC velocities are well calibrated.

Figure 5 shows the dependence of the ${\rm{\Delta }}{v}_{r}$ on position in declination, computed from the maps of Figure 3 weighting, in R.A., on stellar density. The largest offset is found in the central part of the cloud, at $-7^\circ \lesssim \delta \lesssim -6^\circ$, around L1641-N and NGC1980. Here the mean and mode of the velocity distribution have the largest difference, suggesting a very asymmetric distribution. As anticipated and discussed in Paper IV, this part of the Orion A stellar population is $\sim 40 \%$ older than the ONC and L1641-S and is far less extincted, suggesting that it is slightly in the foreground and less strongly associated with the molecular material. A bias due to incompleteness at high A_V is not the case, since the column density of the molecular material in this region is relatively low, up to a few magnitudes in A_V (Lombardi et al. 2014). Thus, star formation in Orion A initiated earlier in the central part of the cloud, which is closer along the line of sight that the remaining gas in the same direction.

Zoom In Zoom Out Reset image size

An offset in v_r between stars and gas may have different origins. One possibility is that the central part of the cloud presented a primordial substructure along the line of sight, both spatially and kinematically. Gas that formed the older population of L1641-N and NGC1980 was closer and with small 1 km s⁻¹ blueshifted motion than the remaining gas today. Such an offset is in line with the typical turbulent motions within a molecular cloud (e.g., Hennebelle & Falgarone 2012). However, a comparison of mean velocities of ¹²CO(2−1) and ¹³CO(2−1) (expected to trace lower and higher density regions, respectively) finds only modest differences of ∼0.1 km s⁻¹. However, in studies of infrared dark cloud (IRDC) filaments, Henshaw et al. (2013) have found offsets of similar magnitude between denser gas traced by N₂H⁺(1−0) and lower-density gas traced by C¹⁸O(1−0).

Another possibility is that prolonged stellar feedback of this foreground population pushed back on the remaining background gas, which is now receding 1 km s⁻¹ with respect to the stellar population.

If primordial substructure along the line of sight is present, the asymmetry of the velocity distribution could also originate from the superposition of two kinematically separated substructures of the cloud. Since, as shown in Paper IV, we found that the v_r distribution around L1641-N does not vary with stellar age, these structures would have needed to sustain star formation at the same time, and with similar duration, which seems somewhat unlikely.

The presence of kinematic substructures in the stellar population of Orion A, and their comparison with those of the gas, can provide important clues on the primordial substructure of the parental star-forming gas. On the other hand, smooth distributions in space and velocity could indicate that dynamical evolution has occurred, since spatial and kinematic substructure is expected to be erased over time in a gravitationally bound system. Evidence of this has been found in the ONC (Da Rio et al. 2014b), where the population in the denser more dynamically evolved core has less spatial substructure than the cluster outskirts.

An analysis of the kinematic substructure in position–position–velocity (ppv) in Orion A was presented by Hacar et al. (2016), using a preliminary version of our IN-SYNC data. These authors adopted a friends-of-friends algorithm that isolates groups or chains of sources in which the separation of each member of the group, in the ppv space, from another member of the group is smaller than a given threshold. Imposing groups with at least 4 sources, they find 37 of such groups, most of them (30) in the northern region of the cloud (L1461N and above), significantly more than expected from random positions in ppv space.

There is no optimal method to identify structures of discrete points in three-dimensional space, rather, different techniques can be developed depending on what one defines as a "structure." In Sections 3.1 and 3.2 we describe two separate methods and discuss the results.

• 1.  
  ppv peaks: this technique is based on building a three-dimensional density map in ppv space from our discrete stellar sample. On this map, then, peaks are identified as overdensities in ppv space. Individual stars can be associated with one or another peak. This method is somewhat analogous to the algorithm Clumpfind (Williams et al. 1994), which is widely used to identify ppv structures in continuous molecular line spectral cube data. This method is sensitive to centrally concentrated overdensities in ppv, but does not necessarily identify filamentary or very elongated structures.
• 2.  
  Connected structures: this technique does not adopt a density map in ppv, but instead builds groups of stars in ppv starting from the highest local density of stars in ppv and adding individual stars as nearest neighbors until a given density threshold is encountered, before creating a new separate group. While this technique does not require kinematic structures to be centrally concentrated, the number and size of structures can be very sensitive to small perturbations of the position of individual sources in ppv space.

3.1. ppv Peaks

To construct stellar density maps in ppv space, we consider each point of the ppv space and find the closest n stars, in ppv, to that point. We adopt a conversion metric such that $1^\circ$ in spatial distance in R.A. and decl. corresponds to 4 km s⁻¹ in v_r. This has been chosen so that the overall spatial extent of the region (i.e., in declination) is comparable with the broadness of the v_r distribution. The stellar density associated with the point is then n/V, where V is the volume of the sphere of radius equal to the distance to the nth star. Therefore, our density maps naturally trace the density of groupings of stars at a given number scale defined by n. Last, coherent structures in ppv space are identified as local density maxima on these density maps within a kernel size of 10' or 0.66 km s⁻¹. For each structure, we also measure the density $\rho (\alpha ,\delta ,{v}_{r})$, as well as the density "contrast" $\delta \rho$ against the local background in ppv space, for which we adopted the density maps obtained for ${n}_{\mathrm{bg}}=30$.

Figure 6 (upper panels) shows the result for n = 5. It is noteworthy that the majority of these groupings are located along the filament south of the ONC, where about half of our IN-SYNC members are concentrated. On the other hand, the ONC as a whole, accounting for almost half of the stellar population in the cloud, is found to be a single kinematic substructure in ppv. Da Rio et al. (2014b) also showed that the projected spatial positions of young members in the ONC, especially the central regions, have reached a relatively smooth quasi-isotropic distribution, and that the population within the half-mass–radius ${r}_{h}\lesssim 1$ pc is at least 7 free-fall times old. This evidence, together with the lack of kinematic substructure we have now identified, points to a relatively old dynamical age for the ONC (Tan et al. 2006), suggesting that primordial structures in space and velocity may have been erased by dynamical interactions. This is in agreement with estimates that the ONC is at least 4 free-fall times older within the half-mass–radius ${r}_{h}\simeq 1$ pc (Da Rio et al. 2014b).

Zoom In Zoom Out Reset image size

Figure 6 also shows that at the location of L1641N and NGC1980 ($-7\buildrel{\circ}\over{.} 5\lesssim \delta \lesssim -6^\circ$) substructure is significant: the different substructures show a large difference in velocity, more than elsewhere. In this region, as shown in Section 2, the overall stellar velocity distribution is asymmetric, thus this asymmetry is likely caused by the superposition of kinematically distinct groups of stars in the same region of the cloud.

After identifying our kinematic substructures as density peaks in ppv, we can attempt to associate individual stars with each substructure. The criterion we chose is the following: a star belongs to a kinematic substructure if, along a segment in ppv space spanning from the star to the density peak, the ppv density increases monotonically. An illustration of the result is shown in Figure 6, bottom panels.

Increasing n decreases the overall number of identified substructures in ppv space, as individual small groupings of stars merge into larger single structures (see Table 1). For n = 20, the only structures detected are three groups for L1641S, two for L1641N (all roughly coinciding with the densest groupings of stars in R.A.–decl. space), plus two structures in the ONC region, and one north of it.

Table 1.  Substructures in ppv Space

	Measured		Random Simulation
n	${N}_{\mathrm{structures}}$	$\delta \rho$	${N}_{\mathrm{structures}}$	$\delta \rho$
4	20	16.347	14.760 ± 2.619	1.405 ± 1.336
5	17	10.740	13.920 ± 3.081	1.184 ± 1.243
6	14	9.215	12.360 ± 2.797	0.937 ± 0.729
7	14	6.486	11.880 ± 2.862	0.899 ± 0.75
8	13	5.427	10.520 ± 1.475	0.851 ± 0.537
10	12	3.431	9.920 ± 1.605	0.702 ± 0.406
20	8	2.050	7.440 ± 1.635	0.500 ± 0.163

Download table as:  ASCIITypeset image

We performed tests to evaluate the statistical significance of our detected kinematic substructures, in particular for what concerns the velocity axis of ppv space. To this end, we created artificial stellar distributions in ppv space, where the positions are the same as our IN-SYNC members, and the velocities are randomly drawn from the local v_r distribution. This latter was computed, for each position in $(\alpha ,\delta )$ measuring the weighted v_r mean of the nearest 30 members, as well as the measured dispersion, after one iteration of sigma clipping, with a $3\sigma$ threshold, to remove outliers. From each realization of artificial distributions in ppv, density maps are created and kinematic substructures isolated, with the same method as was applied on the actual data. For each value of n, 50 of these artificial tests have been performed, and the results are summarized in Table 1. We find that regardless of n, the number of substructures we detect is compatible with the number from random artificial experiments. This, however, is in part the case because our artificial tests kept the stellar positions, and randomized only the velocities. Therefore spatial groupings of sources in the Orion A may still lead to the detection of structures in ppv in spite of a high degree of randomness in v_r. On the other hand, the density contrasts $\delta \rho$ resulting from our data are significantly higher than those from artificial ppv distributions; this confirms that the kinematic substructures we identify in ppv have a statistically significant higher coherence in velocity than the "local" velocity distribution of the nearest 30 stars for each position in the Orion A cloud.

3.2. ppv Connected Structures

The method described in Section 3.1 is particularly sensitive to structures that appear as centrally concentrated "clumps" in ppv space; coherent structures with more irregular shapes may not be detected. Thus, we introduce a conceptually different approach and analyze the different findings. In this case, we do not build a three-dimensional map of stellar density in ppv, but we consider the stellar positions in such space as discrete elements. We associate with each star a local density in ppv as above by measuring the volume containing the closest n stars. We then begin by considering the star with the highest density peak, and use this as the starting point for a group of stars in ppv. We expand this group by adding iteratively one star at a time in order of closest distance in ppv space from any of the group members. The growth of the group halts when the next star about to be added falls at a density below a fixed threshold, ${\rho }_{0}$. Next, a new group is started from the highest density star among those not yet included in a previous group, assuming it is above this threshold.

We vary three parameters: n = 5 or 10, ${\rho }_{0}=1$, or 0.3pc⁻² km⁻¹ s², and also the metric to measure the distance in ppv as the interval ${\delta }_{v}$ on the v_r axis, in km s⁻¹, corresponding a spatial interval of 1°; in this case, we assume ${\delta }_{v}=4$, as in Section 3.1, but also test values of 2 and 8. An example of the results is shown in Figure 7. Overall, we find a higher number of structures than with the method of Section 3.1, as we are able now to identify smaller structures that would not have been resolved given the kernel size in ppv we previously adopted. As a result of this, while the ONC remains the largest identified structure, some additional smaller substructures are now detected in its vicinity. As with the identifications of peaks in the ppv density maps, even in this case, we find the L1641 region to be highly substructured. Decreasing ${\rho }_{0}$ increases the fraction, f, of the total stellar sample included in some identified ppv groups; however, we do not detect a clear difference in the total number of groups. No significant difference in the results is found by changing n or ${\delta }_{v}$.

Zoom In Zoom Out Reset image size

As before, we run the algorithm to artificially generate distributions in ppv space (see Section 3.1), and investigate different diagnostics to test the significance of the results. This is summarized in Table 2. For each parameter choice we compare the number of groups ${N}_{\mathrm{groups}}$ isolated from our data with the mean and standard deviation of this quantity from the random experiments; ${\sigma }_{{N}_{\mathrm{groups}}}$ in this context is the measured excess, in standard deviations, from the result from random experiments. Similarly, ${\sigma }_{f}$ is the departure in standard deviations for the fraction of stars belonging to some structure in ppv. Last, we use a Kolmogorov-Smirnov (KS) test to verify if two cumulative distribution functions are significantly different between our measured values and those from randomly generated distributions. Specifically, we list p₁ as the KS probability (where a low number indicates a significant difference) of the distribution of peak local densities in ppv of the identified groups, and p₂ that of the distribution of number of stars in the detected ${N}_{\mathrm{groups}}$. Results show that for high values of ${\rho }_{0}$, the statistical significance is relatively poor or absent. Decreasing ${\rho }_{0}$, thus enabling a larger number of sources to be labeled as clustered, especially in regions of low density such as L1641, leads to higher significance. Specifically, we identify a number of ppv structures 2–4σ larger than in random experiments, a fraction $f\sim 2\sigma$ higher, and probabilities p₁ and p₂ lower than 1% in most cases.

Table 2.  Groups of Stars Identified in ppv Space and Comparison with Random Tests

${\rho }_{0}$	n	${\delta }_{v}$	${N}_{\mathrm{groups}}$	${\sigma }_{{N}_{\mathrm{groups}}}$	f	${\sigma }_{f}$	p1	p2
1	5	2	36	0.4	30.22	0.87	0.285	0.372
		4	35	0.3	31.2	0.79	0.84	0.105
		8	36	0.6	30.5	0.76	0.291	0.204
	10	2	27	−0.5	38.19	−0.76	0.329	0.626
		4	25	−0.7	41.52	−0.53	0.183	0.003
		8	23	−1.2	45.26	−1.11	0.019	0.533
0.3	5	2	33	2.6	47.0	2.71	0.003	0.071
		4	30	2.2	51.73	1.93	0.173	0.245
		8	36	3.3	43.42	1.94	0.007	0.018
	10	2	29	2.8	52.62	1.88	0.0	0.008
		4	34	4.7	45.12	1.98	0.002	0.001
		8	31	3.2	49.87	1.68	0.0	0.009

Note. See text.

Download table as:  ASCIITypeset image

The fact that we detect statistically significant overdensities in ppv space by adopting two conceptually different methods (Sections 3.1 and 3.2) is highly indicative that kinematic substructure in Orion remains present, and thus may be a tracer of the original turbulent motions of the molecular material in the parental cloud. Future observations (i.e., multi-epoch) will help improve the v_r measurements, especially enabling the identification of binaries. Proper motion measurements from the final GAIA survey will be also useful to study kinematic substructures in Orion. The degree of kinematic substructure that we measure will also be compared to the results from theoretical/numerical models of star cluster formation under different scenarios (Wu et al. 2017a, 2017b).

Limited to the ONC region of the system, we find evidence for a correlation between radial velocities and extinction, in that sources with higher A_V are statistically more redshifted than less deeply embedded members. This is shown in Figure 8, where the offset ${\rm{\Delta }}{v}_{r}$ from the local mean radial velocity (computed as in Section 2 and shown in Figure 3) is plotted against A_V, for all members within a circular aperture of $1^\circ$ centered on the Trapezium region that defines the center of the ONC (e.g., Da Rio et al. 2014b). We also restricted our sample to sources with ${T}_{\mathrm{eff}}\lt 5000$ K, as for higher temperatures the stellar parameters—hence A_V, are more uncertain. A linear fit through the data shows an increase in v_r of about ∼0.1 km s⁻¹/A_V. A randomization test through the permutation of the values along the $y-$ axis of this diagram showed that this correlation is significant at a $\sigma =3.1$ level. The correlation and its significance level both remain unchanged when we vary the aperture diameter between 05 and 2°; the significance reduces at smaller apertures because of low number statistics, and at larger distances outside the ONC.

Zoom In Zoom Out Reset image size

Since A_V is a tracer of the depth of a star along the line of sight—even though the distribution and volume density of the dust is spatially inhomogeneous—this result suggests a global expansion along the line of sight. We exclude that this correlation may be a spurious outcome of correlated systematic effects in the extraction of v_r and A_V since in Cottaar et al. (2015) an opposite trend is observed in IC 348 using APOGEE data and identical data reduction and extraction of parameters as in this work. From the relation ${\rm{\Delta }}{v}_{r}/{\rm{\Delta }}({A}_{V})\sim 0.1$ km s⁻¹/A_V it follows that the central 90% of the population along the line of sight, with $0\lesssim {A}_{V}\lesssim 7.4$, shows an overall v_r spread of ∼0.7 km s⁻¹, or a net maximum expansion from the cluster midplane of ∼0.35 km s⁻¹. The reddening distribution is highly asymmetric, with a median A_V—which could be associated with the cluster midplane—of 1.3 mag. Note that this fact is not only seen in our IN-SYNC survey, but in all previous studies reaching higher completeness and A_V (Robberto et al. 2010; Da Rio et al. 2014b, e.g.,); thus it is not the result of completeness limits of this work. This implies that the volume density of the molecular material increases along the line of sight. We are unable to explicitly measure the net velocity offset of two halves of the population from this midplane, as in the low-extinction limit, we do not measure any correlation of v_r with A_V. This, however, is likely because the typical uncertainty in A_V is ∼0.4 at all extinctions. Thus, for low A_V , the relative uncertainty is particularly high, flattening any v_r–A_V correlation due to regression dilution.

We look for variations in the measured correlation indicating expansion with mass and age, and divide the stellar samples into two mass ranges ($M\lt 0.3{M}_{\odot }$ and $0.3{M}_{\odot }\lt M\lt 1{M}_{\odot }$) and two age ranges on both sides of the median cluster ages. These parameters have been obtained from the Siess et al. (2000) evolutionary models (IV). Although there is a hint that lower-mass stars show a large velocity gradient with A_V, and so do older stars compared to younger ones (see Figure 9), these differences are not statistically significant ($\sim 1\sigma$). Future data are needed to confirm or refute a mass and age dependence of the expansion of the ONC.

Zoom In Zoom Out Reset image size

Finally, no correlation between radial velocity and extinction is found along the filament south of the ONC (L1641-N and L1641-S).

We obtained the stellar velocity dispersion ${\sigma }_{v}$ for different samples of the Orion A population in a similar fashion as in Paper II and III for NGC 1333 and IC 348, respectively, in Perseus. Specifically, we consider the observed distribution of v_r as well as the measurement uncertainties in radial velocity, and the contribution from unresolved binarity. Unlike the IN-SYNC survey in Perseus, our v_r values have only one epoch for the majority of our stellar sample, therefore our binarity correction has been performed statistically. We adopt a flat mass ratio distribution (Reggiani & Meyer 2013), a flat eccentricity distribution (Duchêne & Kraus 2013) over the whole mass range, and the log-normal period distribution from Raghavan et al. (2010). Using the tool VELBIN (Paper II), we then generated a large number of binary orbits randomly drawn from these parameter distributions, and assuming random orientations, we obtained the one-dimensional distribution of radial velocities of the primary. We assumed a binary fraction ${f}_{\mathrm{bin}}=0.44$, as determined by Raghavan et al. (2010) for solar-type stars.

Considering a sample of v_r and associated errors $\delta {v}_{r}$, we generate a large number of intrinsic velocity dispersions ${\sigma }_{v}$ as normal distribution of varying standard deviation. We then convolve these with the contribution due to measurement errors; to this end, we sum the individual distributions of errors, all assumed normal, for each star. Since stars have different errors, the resulting distribution is not a Gaussian. Last, we consider the one-dimensional distribution of v_r offsets from binaries for ${f}_{\mathrm{bin}}=1$, scale it down to our value of ${f}_{\mathrm{bin}}$, complement it with $(1-{f}_{\mathrm{bin}})$ zeros representing single stars, and convolve it to the intrinsic ${\sigma }_{v}$ plus measurement error distributions. The resulting simulated distributions as well as the measured v_r distribution are fitted with Gaussian distributions after removing outliers at $3\sigma$, and the best match is found. This naturally represents the best velocity dispersion ${\sigma }_{v}$, which, convolved with contributions from errors and binaries, reproduced the observed v_r distribution. An example of this procedure is shown in Figure 10 for a sample of members centered on the ONC within a radius of $0\buildrel{\circ}\over{.} 3$ or ∼2.2 pc.

Zoom In Zoom Out Reset image size

Generally, the contribution of the measurement errors is far smaller than both the observed v_r distribution and the intrinsic ${\sigma }_{v}$. For example, the distribution of v_r offsets due to measurement errors shown in Figure 10 is such that 68% of the data (i.e., equivalent to a 1σ interval) is in the range ±0.38 km s⁻¹. Furthermore, as pointed out in Paper III, the contribution from binaries is even smaller (68% of v_r offsets in the range ±0.10 km s⁻¹ for ${f}_{\mathrm{bin}}=0.44$). Therefore our results are not very sensitive to the assumptions made for the binary population, i.e., the distribution of binary parameters and ${f}_{\mathrm{bin}}$. This fact would be less the case for high-mass stars, both because of a higher binary fraction, and possibly because of a lower intrinsic ${\sigma }_{v}$ from dynamical energy equipartition. However, our sample is dominated by low-mass stars, with 95% of members with a mass lower than 2 M_⊙.

Figure 11 shows the derived corrected ${\sigma }_{v}$ in different bins in declination along the Orion A cloud. The solid line represents the dispersion as is, adopting the measured v_r for all stars in each bin. We have also corrected for the spatial gradients in v_r along the region (see Section 2), which if not accounted for, may lead to an overestimate of the overall ${\sigma }_{v}$ in a region. To this end, we adopted the radial velocities as the difference, star by star, between the measured v_r and the local mean radial velocity in that position, as computed in Section 2 and shown in Figure 3, second panel. We refer to these two dispersions as "standard" and "gradient-corrected." The result is shown in Figure 11 as a dashed line. Overall, ${\sigma }_{v}$ varies in the range 1–2.5 km s⁻¹. The most striking feature is a broad peak in ${\sigma }_{v}$ in L1641N; this is consistent with the detected asymmetry of the velocity distribution and the presence of kinematic substructure in this region, possibly indicating spatial structure along the line of sight (see Section 3). Thus a larger ${\sigma }_{v}$ could be the result of a superposition of stellar systems with relative bulk mean motions. Just south of the ONC, we detect the most prominent drop in velocity dispersion. The ONC also reveals itself as a local peak in ${\sigma }_{v}$ at values 2–2.5 km s⁻¹. Moving north of the ONC, in the region of NGC 1980 and NGC 1977/OMC 2–3, ${\sigma }_{v}$ remains roughly constant at values ∼2 km s⁻¹. These results are in fair agreement with the optical radial velocity survey of Tobin et al. (2009), which covered a declination range $-6^\circ \lt \delta \lt -4^\circ$. The velocity dispersion of ¹³CO gas is also shown in Figure 11; the stellar ${\sigma }_{v}$ is higher than that of the gas the regions of L1641N and around the ONC, whereas between these two, as well as in the southern end of the filament (L1641S), stars show velocity dispersions that are in more agreement with the CO.

Zoom In Zoom Out Reset image size

Next we focus in more detail on the ONC. Figure 12 shows the radial dependence of ${\sigma }_{v}$ from the center of the cluster. Again, we compute both the "standard" and "gradient-corrected" dispersions. These average at ∼2.1 and ∼1.7 km s⁻¹, respectively. We note, however, that given the small spatial size of the ONC compared to the Orion A cloud, or the width of the declination bins in Figure 11 (3 pc correspond to ∼04), large-scale gradients in v_r have a smaller effect on the measured velocity distribution. Still, we do compute the "gradient-corrected" ${\sigma }_{v}$ along with the "standard one." Our measured standard ${\sigma }_{v}$ is in fairly good agreement with the dispersion of ∼2.3 km s⁻¹ measured in the ONC by of Tobin et al. (2009), as well as the stellar proper motion study of Jones & Walker (1988); none of these studies, however, had derived ${\sigma }_{v}$ after accounting for spatial variations in the mean ${\sigma }_{v}$, hence our lower value for the gradient-corrected dispersion.

Zoom In Zoom Out Reset image size

In Figure 12 we also show the predicted velocity dispersion profile required for dynamical virial equilibrium, as derived in Da Rio et al. (2014b) from the total mass profile of the ONC. This was obtained by combining estimates of the total stellar mass profile from censuses of stellar parameters from the literature, complemented by X-ray and infrared surveys, and estimates of the gas volume density within the ONC. Specifically, the stellar component describes a power-law profile slightly steeper than a single isothermal sphere, and the gas was assumed to be at a constant volume density of 22 M_⊙ pc⁻³. This gas density was derived from the stellar A_V distribution, which averages at ${A}_{V}\sim 3$ mag, and not from, e.g., the total column density of molecular material; this is because the vast majority of the ISM along the line of sight of the ONC, reaching up to 100 mag (Bergin et al. 1996; Lombardi et al. 2014), lies behind the cluster. The gas contribution is negligible compared to the stellar mass up to a few half-mass radii (${r}_{h}\simeq 1$ pc, Da Rio et al. 2014b), and is responsible for the increase of the equilibrium ${\sigma }_{v}$ at increasing distances from the center. The mean ${\sigma }_{v}$ in the ONC required for virial equilibrium based on the total mass is ∼1.73 km s⁻¹ (Da Rio et al. 2014b), which, compared to the measured 2.2 km s⁻¹, indicates that the ONC is supervirial with a virial parameter (Bertoldi & McKee 1992) ${\alpha }_{\mathrm{vir}}\simeq 1.5$. On the other hand, our gradient-corrected ${\sigma }_{v}$ indicates a fully virialized system. The predicted virial velocity is not without uncertainties: in Figure 12 we also show the range in the radial profile of the expected virial velocities assuming an uncertainty of 30% in the estimated mass. This uncertainty is likely quite large for the stellar component, being derived from deep samples at different wavelengths. On the other hand, a higher uncertainty may come from the gas component, especially if the distribution of the material along the line of sight—hence the conversion from column densities to volume density—is quite different than in the orthogonal directions, specifically, if the system is more flattened along the line of sight.

Our results show that the ONC population appears in very good agreement with a virial state when considering the gradient-corrected ${\sigma }_{v}$, i.e., correcting for local variations in the mean v_r. On the other hand, adopting the standard ${\sigma }_{v}$, stars are slightly supervirial. Even in this case, however, velocities are well into the boundness regime (${\sigma }_{\mathrm{bound}}=\sqrt{2}{\sigma }_{\mathrm{vir}}$).

Even though our results are consisten with a virial state for the ONC, we comment on the scenario of a moderately supervirial scenario. This would be compatible with our presented evidence that the cluster is currently expanding (see Section 4); however, a small population of expanding stars can be present even with the bulk of the cluster being virialized. This could be attributed to some missing mass from a recent gas expulsion from the cluster location. Numerical simulations (e.g., Farias et al. 2015), however, have shown that the virial state of a cluster fluctuates significantly during the early evolution of a stellar system, and the fraction of the stellar population that remains subsequently bound strongly depends on the dynamical state at the time of gas removal. Additional data, such as a scheduled continuation of the IN-SYNC Orion survey, to increase the stellar sample and obtain additional epochs for our IN-SYNC targets, as well proper motions from HST photometry (HST program GO13826, PI M. Robberto) or GAIA for the bright end of the population, will allow us to better characterize this potential expansion of the ONC. It is also important to stress that a possible infall of material (gas and stars) along the filament may affect the future fate of the system. In fact, Tobin et al. (2009) suggested that the large velocity gradient, for both stars and gas, that is observed to the north of the ONC (see Figure 1) could indicate an infall of material toward the ONC.

On the other hand, the ONC being close to or fully virialized, is compatible with the evidence that the system is dynamically evolved (Da Rio et al. 2014b). The fact that ${\sigma }_{v}$ for stars in the ONC is larger than the molecular gas dispersion (see Figure 11) is not of concern. This is because, as anticipated earlier, most of the molecular material in this region of the cloud is located behind the young stellar population, as confirmed by the very low extinction affecting the cluster. We do not find variations of ${\sigma }_{v}$ as a function of stellar mass in the ONC. However, we stress that the IN-SYNC survey does not cover massive stars, as the sample is limited to below a few ${M}_{\odot }$. We do measure a slight variation with age: stars younger than the median age show a ${\sigma }_{v}$ 0.2 km s⁻¹ smaller than older stars; this is significant only at a 1.5σ level, however, and could originate from biases. On the one hand, older PMS stars are fainter and thus affected by larger errors; small systematics in the estimate of the errors in v_r may therefore be responsible for the detected difference. Moreover, as shown in Da Rio et al. (2012), at fainter luminosities contamination from background sources increases. Even though we restricted our analysis to known or probable members, these membership estimates are heterogenous and based on methods with a range of reliabilities in excluding false candidates.

Returning to the L1641, south of the ONC, Figure 11 showed that the stellar velocity dispersion in the northern region (L1641N) is in the range 2–2.5 km s⁻¹, not very dissimilar from that of the ONC, whereas ${\sigma }_{v}$ drops to values ∼1 km s⁻¹ in L1641S. In both these regions, the total density is much lower than that in the ONC, so that the stellar population should not be bound by its own mass. The average stellar density at declinations $\delta \lt -6^\circ$ is ∼10 pc⁻². When we assume a depth of the population along the line of sight comparable to the width of the filament in the plane of the sky ($\sim 1\buildrel{\circ}\over{.} 5$ or 10 pc), this results in a stellar number volume density of $n\lesssim 1$ pc⁻³ or ${\rho }_{* }\lt 1\ {M}_{\odot }$ pc⁻³. This is more than two orders of magnitude lower than in the ONC, where ${\rho }_{* }\simeq 250\ {M}_{\odot }$ pc⁻³ within r_h.

As for the gas density, we adopt the extinction map from Lombardi et al. (2014); within the Orion A southern filament ($\delta \lt -6^\circ$), the mean extinction is ${A}_{V}\simeq 1$ mag. As shown in Da Rio et al. (2014b), this corresponds to a mass column density of $\sim 22\ {M}_{\odot }$ pc, which again, assuming a depth of the filament of ∼10 pc, results in a volume density of molecular material $\sim 2\ {M}_{\odot }$ pc⁻³.

We compare these estimates with the condition for virial equilibrium of a filamentary structure from Fiege & Pudritz (2000), i.e.,

$$\begin{eqnarray}&&{\sigma }_{{vir}}\sim \sqrt{\displaystyle \frac{{{Gm}}_{l}}{2}},\end{eqnarray} \tag{ 1 }$$

where m_l is the mass per unit of length. Considering a width of the L1641 filament of $\lesssim 10$ pc, we obtain a typical ${m}_{l}\sim 200\ {M}_{\odot }$ pc⁻¹, resulting in ${\sigma }_{\mathrm{vir}}\simeq 0.65$ km s⁻¹. This is smaller than our measured ${\sigma }_{v}$ throughout the L1641 region. We perform a more rigorous analysis of these quantities in different positions of the filament. This is shown in Figure 13, where the displayed quantities, as a function of declination, have been computed within moving strips of 1 pc in width in the plane of the sky perpendicular to the filament. The width of the filament has been measured by imposing lateral boundaries where the extinction drops to A_V = 0.3 mag according to the map of Lombardi et al. (2014). For the stellar masses at each position, we adopted the values derived in Paper IV; for guidance, we also include the results from doubling the stellar masses, as an extreme case to account for residual incompleteness in our IN-SYNC sample. This is largely irrelevant, as the total mass is gas dominated.

Zoom In Zoom Out Reset image size

The value ${\sigma }_{\mathrm{vir}}\sim 0.6$ km s⁻¹ is roughly constant along the filament; such a virial dispersion is several times smaller than the observed ${\sigma }_{v}$ in L1641N, indicating that this region is highly unbound. As shown in Paper IV, here the dust extinction affecting stars is ${A}_{V}\lt 1$ mag, and the velocity dispersion is highly asymmetric (Section 2). Therefore the population here is mostly in the foreground of the remaining molecular material, and possibly populated by substructured groups along the line of sight, showing an overall velocity dispersion consistent with being unbound. On the other hand, ${\sigma }_{v}$ in L1641S, while being somewhat higher than the virial one, is not incompatible with it, given the uncertainties.

Assuming the present ${\sigma }_{v}$, the crossing time ${t}_{\mathrm{cross}}$ (see Figure 13) ranges from ∼3 Myr in L1641N to ∼6 Myr in L1641S. This is somewhat older than than the age of the system and is therefore compatible with the fact that the young population remains fairly well spatially associated with the molecular material. As we emphasized above, however, our poor knowledge of the distribution of material, and in particular molecular gas, along the line of sight may affect our estimate for the virial equilibrium velocity predicted for L1641.

If the filament is regarded as being a single dynamic entity, then this can explain why the ¹³CO velocity dispersion is quite close to the expected viral velocity dispersion. However, then it unclear how the stellar component would have gained such a higher kinetic energy. First, the relatively high dynamical time in the region grants that the stellar population is dynamically young, thus stellar interactions and binary star evolution are nor responsible for the observed velocity dispersion. Also, the L1641 region is deficient of massive stars (Hsu et al. 2012, 2013), thus a rapid removal of gas from stellar feedback should be ruled out. Stutz & Gould (2016) recently proposed a model that aims to explain how velocities of young stars can be larger than that of dense gas. They noted that young, but already-formed, stars in Orion A appear more spread out than forming protostars both in the plane of the sky, as well as in their radial velocities (based on IN-SYNC data). The protostars are spatially and kinematically associated to the filaments of molecular gas. They suggested that magnetically induced tangential oscillation of the filament may be taking place. In this scenario, stars are formed coupled with the dense gas, and then gently ejected as the molecular material oscillates back to its original position. However, it is difficult to see how this scenario could lead to observed line-of-sight velocity dispersions of the ¹³CO traced dense gas that are consistently much smaller than those of the ejected stars along the entire filament.

We have continued our analysis of the IN-SYNC Orion survey, a high-resolution near-infrared spectroscopic survey of about 2700 stars in the Orion A star-forming region, that was obtained with the SDSS APOGEE spectrograph. It covers $\sim 6^\circ$ in the plane of the sky, or ∼40 pc, thus covering the entire cloud from NGC 1977, the ONC, L1641N down to L1641S. We cover a mass range from $\sim 0.1\ {M}_{\odot }$ to a few solar masses and extinctions up to about 20 mag. In the previous paper of this series (IV), we presented the derived stellar parameters, studied the age and age spread of the system, the reddening law, and assigned new membership estimates for sources without a previous membership confirmation.

Here we have focussed on the stellar kinematics and dynamics in the region, from the analysis of radial velocities. Here we summarize the main results.

• 1.  
  We find that young stars have average radial velocities that are very similar to those of the molecular gas throughout the region as traced by CO, following a ∼10 km s⁻¹ variation from north to south, with the northern region redshifted and the southern region blueshifted.
• 2.  
  In the central part of the cloud, corresponding to L1641N, stars appear slightly blueshifted with respect to the gas. In this region, however, stars are older and affected by little extinction, suggesting that the molecular gas is located behind the population along the line of sight, and receding with respect to the stars.
• 3.  
  In the region of L1641N, the distribution of radial velocities v_r appears asymmetric, with a broader blueshifted tail. We interpret this as a possible superposition of stellar substructures along the line of sight in relative motion with each other.
• 4.  
  We find spatial and kinematic substructures in the young population throughout the region, as statistically significant overdensities in position–position–velocity space. These may indicate groupings of stars that have preserved the turbulent motions of their natal gas. In the L1641N region, these structures show the largest relative offsets in radial velocity, corroborating the possibility of a highly structured stellar distribution along the line of sight. The ONC, on the other hand, does not present significant kinematic substructure, in agreement with previous findings for this system, which reported that it has undergone some degree of dynamical evolution, which is expected to lead to phase mixing of stellar orbits.
• 5.  
  We find indications that the ONC may currently be expanding, from the correlation between radial velocities and A_V, with the low-extinction end of the population, which is closer to the observer along the line of sight, blueshifted, and, on the contrary, the high-extinction stars redshifted.
• 6.  
  In the ONC, we measure a velocity dispersion ${\sigma }_{v}\simeq 2.2$ km s⁻¹, when considering the measured v_r, in agreement with previous findings in the region based on radial velocities or proper motions. On the other hand, when correcting for spatial variations in the mean v_r—which is critical for large-scale gradients in v_r, but does marginally affect the results at small scales—${\sigma }_{v}$ decreases to ∼1.7 km s⁻¹. Comparing these values with estimates of the total mass in the system, they result in the ONC being slightly supervirial (but bound), or in good agreement with a virialized state, respectively. Therefore we cannot exclude that the ONC has reached virial equilibrium.
• 7.  
  The southern part of the Orion A cloud, L1641, shows two separate behaviors. The northern part, L1641N, has a velocity dispersion in the range 2–2.5 km s⁻¹, similar to that of the ONC. However, considering that the density of stars and gas is much lower (a few ${M}_{\odot }$ pc⁻³) than in the ONC, this is a few times higher than the velocity dispersion required for (filamentary) virial equilibrium (∼0.6 km s⁻¹), pointing to an unbound population. The low A_V in this region suggests that stars have decoupled from the gas, which is located mostly in the background and had velocities consistent with a virial. In the southern end of the region, L1641S, on the other hand, ${\sigma }_{v}$ is lower and in approximate agreement with (or slightly greater than) both a virial state and the velocity dispersion of the gas.

J.J.T. acknowledges support from the University of Oklahoma, the Homer L. Dodge endowed chair, and grant 639.041.439 from the Netherlands Organisation for Scientific Research (NWO).