\textcolorblackGravity or turbulence V: Star forming regions undergoing violent relaxation.

In the following sections we analyse six different numerical simulations that represent the evolution of a molecular cloud. The first two where performed with GADGET-2 (Springel, 2005), which is a smoothed particle hydrodynamics (SPH) code, with the inclusion of sink particles as proposed by Jappsen et al. (2005). The remaining four were performed using Flash 4.3 (Fryxell et al., 2000), with the inclusion of sink particles as proposed by Federrath et al. (2010). These simulations represent different physical conditions of molecular clouds, according to the current models of cloud formation and evolution. We briefly describe these simulations, quoting the corresponding papers for details of each run.

We have derived stellar masses and ages for a sample of kinematic candidates selected using the astrometric observables of parallax ($\varpi$) and proper motions ($\mu_{\alpha}$, $\mu_{\delta}$) from GAIA-EDR3 (Gaia Collaboration et al., 2021) in the chosen studied region (hereafter, the ONC region: $83.0^{\circ}\leq\alpha_{J2000}\leq 84.5~{}\circ$ and $-6.5^{\circ}\leq\delta_{J2000}\leq-4.0^{\circ}$). We apply the selection criteria used by Hernandez et. al. (in preparation), which we describe here for clarity:

1. 1.

   We first apply the zero-point correction of parallax following the method described by Lindegren et al. (2021b).

2. 2.

   We restrict the astrometric quality of the sample, requiring a Re-normalised Unit Weight Error (RUWE)$<$1.4 (Lindegren et al., 2021a) and a parallax error smaller than 20%.

3. 3.

   We then selected stars with parallaxes in the range $2.27<\varpi<2.93\;\rm mas$ and stars with a proper motion modulus ($\mu=\sqrt{\mu_{\alpha}^{2}+\mu_{\delta}^{2}}$) $<3.23\;\rm mas\,yr^{-1}$. Based on the $\varpi$ and $\mu$ distributions of the kinematic young members selected by Kounkel et al. (2018a) in the the ONC region, these limits were defined using the median and the standard deviation applying a 3$\sigma$ criteria.

With these filters, we obtain a sample of 3030 kinematic candidates. We first locate these stars with high precision in the H-R diagram, and use models of pre-main-sequence stellar evolution, in order to infer their masses and ages. For this purpose, we cross-match our kinematic candidates with the optical spectroscopic census of Hernandez et al. (hereafter the optical spectroscopic sample; in preparation), finding spectral types for 1586 stars. Then, we used the MassAge code (Hernandez et al., in preparation) to compute ages and masses of 1461 stars. The sample without masses and ages includes the star 2MASS_J05343988-0625140, which falls below the main sequence, and 124 stars which do not have reliable $J-$band magnitude.

In brief, the MassAge code uses spectral types (or effective temperatures), photometry ($G_{p}$, $R_{p}$, and $B_{p}$) and parallaxes from GAIA-EDR3 (Gaia Collaboration et al., 2021), and the $J$ and $H$ magnitudes from 2MASS (Cutri et al., 2003). The uncertainties in the estimated values are obtained using the Monte Carlo method of error propagation (Anderson, 1976). The extinction is obtained by comparing the observed colors with the standard colors reported in Luhman & Esplin (2020). Here, we use the extinction law from Cardelli et al. (1989) to obtain the extinction in each photometric band normalized at 0.55µm ($A_{\lambda}/A_{V}$). The luminosity is estimated from the $J-$band absolute magnitude using the bolometric correction from Pecaut & Mamajek (2013). We convert spectral types into effective temperatures using the standard table of Pecaut & Mamajek (2013). Finally, we obtain stellar masses and ages by comparing the location of the stars on the Hertzsprung–Russell diagram with theoretical values. Here, we use two evolutionary models: the PARSEC model (Marigo et al., 2017) and the MIST model (Dotter, 2016).

In order to expand our sample, we have included 382 and 146 additional stars studied by Da Rio et al. (2012b) and Olney et al. (2020a), respectively. Da Rio et al. (2012b) reported effective temperatures using two narrow-band filters located at 7530 Å and 7700 Å (tracing the continuum and a TiO-band feature, respectively) and a broad I-band filter, taking into account the extinction law from Cardelli et al. (1989). Additionally, using a deep convolutional neural network, the APOGEE-net tool (Olney et al., 2020a) improves the effective temperature and the stellar surface gravity derived by Kounkel et al. (2018a) by comparing APOGEE H-band spectra (R $\sim$22500) and PHOENIX spectral theoretical library (Husser et al., 2013). On the other hand, the optical spectroscopic sample includes 429 and 697 stars studied by Da Rio et al. (2012b) and Olney et al. (2020a), respectively. Generally, the effective temperatures derived in those works agree within 500K with the effective temperatures derived in the optical spectroscopic sample. Using the effective temperature as input parameter in the MassAge code, we derived stellar masses and ages for stars studied by Da Rio et al. (2012b) and Olney et al. (2020a), not included in the optical spectroscopic census of Hernandez et al. (in preparation). Figure 2 indicates that the additional stars increase the number of faint kinematic candidates ($G>12$) with estimated masses (likely low mass stars). We include a total of 1989 stars with ages and masses (filled histogram), representing more than 65% of the kinematic candidates.

These sources are spread over an area of 7$\times$7 pc, which we call the extended area, or the total sample. In addition, we distinguished the central core, a region of about 1.5$\times$1.5 pc, determined as the overdensity having more than 30 stars per square degree, 3 times the rms value of the extended area. Fig. 3 shows a map of the analysed region. The stars in the extended region are shown as blue circles, and the stars in the central cluster as magenta circles.

Fig. 4 shows the age histograms of our two samples of the ONC. The top panel shows the ages inferred using MIST models, while the lower panels show the corresponding distributions using PARSEC models. Some points are worth stressing. First of all, that each model produces a statistically different age distribution from the other model. For the total sample (blue histograms), for instance, the MIST models produce younger distributions than the PARSEC models, while for the the central sample (red histograms), the situation is reversed and the MIST models are the ones that produce older distributions.

Second, the number of stars older than 5 Myr strongly decreases in the case of the MIST models. In comparison, the PARSEC models have still important contribution of stars for bins up to 10 Myr. We will discuss this and the previous point in §4.2.

Third, we notice that the central sample (red histograms) is statistically younger than the extended or total sample (blue histograms) regardless the evolutionary model, a fact known to occur also in other clusters (e.g., Getman et al., 2014; Getman et al., 2018).

Fourth, in the present work we will consider only the masses and proper motions of the stars younger than 10 Myr in each sample, in order to minimize contamination by older dispersed populations and/or possible stars with edge-on disks that appear older because of disk obscuration. We notice that the contribution of older than 10 Myr stars to the statistics shown in the next sections will be negligible: in the case of the total sample, we are rejecting only 5% (MIST) and 13.7% (PARSEC) of the stars, while in the case of the central core, 2% (MIST) and 9% (PARSEC).

In Fig. 5 we now show the resulting mass histograms of the extended (blue) and central (red histograms) obtained with the MIST (top) and PARSEC (bottom) models. Although the distributions from both models are again statistically different, we will show in §3.4 that our result based on the velocity dispersion of the stars per mass bin is robust, regardless the model used.

Fig. 6 shows the velocity dispersion of the sinks in different mass bins, at $\sim$1$t_{\rm ff}$ for the clusters in our simulations of global, hierarchical and chaotic collapse (runs $1-3$, see §2). As it can be seen, despite the variety of initial conditions in our suite of simulations, all of them exhibit a fairly flat velocity dispersion as a function of the mass bin of the stars. This behaviour is the expected result for a cloud undergoing a dynamical collapse, and thus, suffering a violent relaxation process (Lynden-Bell, 1967).

In contrast, the case of sinks in simulations of clouds with forced turbulence is substantially different. In Fig. 7 we plot the velocity dispersion per mass bin of stars formed after stirring our turbulent box during 5 dynamical crossing times, and right after gravity is turned-on and allows the formation of sink particles. This figure shows that the velocity dispersion increases with the mass of the bin by a factor of $\sim$2 when the mass increases by a factor of 10, suggesting a substantially different process than collapse, which produces ‘hotter’ (in the kinetic sense, i.e., with larger velocity dispersion) massive stars than low-mass stars.

Thus, we can argue that, while collapsing clouds will produce stars that exhibit constant velocity dispersion per mass bin due to the violent relaxation, turbulent environments will produce stars that are kinetically segregated in mass, i.e., with the more massive stars having larger velocity dispersion than the less massive stars.

In order to characterise furthermore the kinematics during the formation of the clusters in our simulations, we define the expansion factor as

where $\hat{r}_{*}$ is the unit vector of the star measured from the center of mass of the cluster, $\mathbf{v}_{*}$ its velocity vector, and $i=1,2,3$ are the Cartesian coordinates, such that the Einstein convention for repeated indexes is adopted. The summation runs over the $N$ stars in the sample.

In Fig. 8 we plot the expansion factor ($\cal{EF}$) as a function of time, for the cluster formed in Kuznetsova et al. (2015, a detailed study of the expansion factors in the full set of simulations is presented in Bonilla-Barroso et al., in preparation). The run was stoped at $t\sim$ 0.7 Myr, one free-fall time of the gas at its initial density $\rho_{0}$ in the box,

(where $G=6.67\times 10^{-8}$ gr cm^-3 s^-2 is the constant of gravity, $n=\rho/\mu m_{H}$ is the number density, $m_{H}$ is the mass of the atom of hydrogen, and $\mu=2.36$ the mean molecular weight, assuming solar abundances in MCs). The dashed line in this figure, labeled as 3D, denotes the exact calculation as defined in eq. (3). However, since proper motions only measure a particular projection of the actual 3D data, observations can be skewed by the particular point of view between the earth and the cluster. Thus, we also computed the 2D version of eq. (3) through 1,000 different random directions, and used only the 2D motions projected in a plane perpendicular to each line of sight. As a result, in Fig. 8 the solid line, labeled as 2D, shows the mean value of the 2D version of the expansion factor (3), while the dark and light shaded areas represent the rms and 3 rms values around the mean, respectively.

Several lessons can be learned from this figure. First of all, it shows us that, at the beginning, the expansion factor appears to be negative, as expected for any collapse process. However, as the potential well increases and more gas and stars fall into the center, the clusters suffer a series of expansions and contractions (we avoid to call them oscillations because that may give the false impression of a stationary system, which is not the case, since the cluster continues accreting gas and stars, and forming additional stars). This is the period of time during which the velocities are seriously randomized (Fig. 7 by Kuznetsova et al., 2015).

Naively, one may think that this randomization can give rise to a more dynamically relaxed configuration. However, the relaxation time $t_{\rm relax}$, defined by eq. (2), is still large compared to the evolutionary time. In the present case, the relaxation time (2) is of the order of $t_{\rm relax}\sim$12$t_{\rm cross}$, since $N_{*}\sim 800$, and increasing with time as more stars are born. Thus, assuming that the time span in each oscillation of the cluster is precisely one crossing time $t_{\rm cross}$, the cluster has not had time to relax. Thus, even though the velocities are randomized, the velocity dispersion of the stars in a cluster formed in a global collapse process is independent of the mass of the stars (Lynden-Bell, 1967), as shown in Fig. 6.

A second and important point to mention is that, statistically speaking, even though the cluster has formed during a global collapse of its parental cloud, the detailed value of the ${\cal EF}$ depends on the particular 2D projections in which it is observed, as the shaded areas in Fig. 8 shows, as well as the exact time in which the core is observed. Thus, attributing a particular dynamical state to a cluster that it is still in the process of formation and has substantial amounts of gas in its surroundings, based on the current proper motions (e.g., Rivera et al., 2015; Dzib et al., 2018; Ward & Kruijssen, 2018; Ward et al., 2020), should be taken with caution.

In Fig. 9 we show the velocity dispersion per mass bin of our sample of stars in the ONC. Left panels contain stars within the $\sim$1.5 pc central ONC, while right panels contains the stars of the more extended ($\sim$7 pc along north-south) region. The masses in the upper panels were computed using the MIST models, while in the lower panels the PARSEC models were used. It is clear from this figure that the velocity dispersion of the stars per mass bin has no dependency with the mass of the bin, a result that spans a factor of $\gtrsim$ 60 in mass. This result appears to be consistent with the scenario of collapse, and clearly calls into question the somehow spread idea that the ONC formed out of a dense core supported by turbulence against global collapse along $\sim$10 free-fall times at the current density, and forming stars during this time.

An important question to address is in order. If the MIST and PARSEC models produce statistically different mass distributions (as shown in Fig. 5), could the none-dependency with mass of the velocity dispersion (Fig. 9) be the result of actual masses of the stars mixed up in different bins, such that actual differences in the velocity dispersion of the stars with different masses is averaged out? The answer is no. As shown in Fig. 10, it is clear that the massive stars in both models are not mixed with the low mass stars, and thus, the more massive bins (above $\sim$1 $M_{\odot}$) contain the same stars in both models, and thus, the same statistics. Any increase or decrease of the velocity dispersion, at least of the most massive stars, will be detected, regardless the model used.

A frequent misconception is that the proper motions of the stars in collapsing clouds are necessarily inward motions, if the observation occurs during the process of contraction of the cloud, and outwards, once the remaining gas cloud has been removed by stellar feedback and the potential well disappears (e.g., Ward & Kruijssen, 2018; Ward et al., 2020). In reality, the hierarchical and chaotic collapse that occurs in molecular clouds is never as simple as such picture. In fact, as collapse proceeds, denser regions keep incorporating more gas and stars from the surroundings, increasing the potential well, and allowing the stars in the cluster to rapidly randomise their velocity vectors as the potential well increases as the result of the infall of material (see Figs. 7 and 13 from Kuznetsova et al., 2015).

Several authors have argued that the ONC should be clearly expanding if it was formed by the GHCC model, because the gas has already been removed by the stars (e.g., Krumholz et al., 2019). We want to stress several points in this regard. First of all, there is a considerable amount of molecular gas that remains behind the optically visible cluster (e.g., O’Dell, 2001). As pointed out by Hillenbrand & Hartmann (1998), a reasonable estimate is $\sim$ 4000 $M_{\odot}$ of gas within a diameter of 4 pc (2200 $M_{\odot}$within 2 pc, see Bally et al., 1987).

In addition, as we have shown in §3.3, during a substantial fraction of the formation timescale of a still embedded cluster, the expansion factor ${\cal EF}$ (eq. [3]) can oscillate around zero, even if the cluster has been formed by the global collapse of its parental clump³³3We recall also that the detailed value of ${\cal EF}$ could depend upon the exact time in which it is observed, and its 2D calculation from proper motions, furthermore, may depend on the detailed projection it has, as discussed in §3.2.. Furthermore, it is well known that there is a population of embedded stars in the region, and thus, the cluster can still be under construction, with stars and gas falling in, as it occurs in numerical simulations (e.g., Kuznetsova et al., 2015).

Finally, it is important to recall that the original problem posed by Krumholz et al. (2019) is not exclusive of the GHCC model. In other words, the fact that ONC exhibits a velocity dispersion that exceeds its virial value (and thus, must be expanding), but has a small expansion factor, is a problem that can be posed to any model of molecular cloud evolution and cluster formation, regardless of whether the original cloud is collapsing, or supported by turbulence.

We conclude, thus, that neither the apparent randomization of velocities in the ONC, or its small rate of expansion, can be an argument to dismiss the global collapse of a larger clump as progenitor of the ONC. The fact that the ONC region still has substantial amounts of gas it its background should be considered when arguing whether expansion should be observed or not.

Its been argued that the stars in the ONC exhibit an age histogram with a peak at $\sim$6–8 Myr (Krumholz et al., 2019, see panels c and d in their Fig. 13), for an area similar to our extended region, and nearly constant for the central cluster. These authors argue that the ONC has been forming stars over the last 6 $-$ 8 Myr, and thus, over $\gtrsim$10 $\tau_{\rm ff}$. Several points need to be addressed in this respect.

There are several similar definitions of the efficiency of star formation, but the most simple one is the amount of gas of a cloud that has been converted into stars. In clusters where stellar feedback has cleared up at least partially the mass of the parent cloud, it becomes necessary to estimate the total mass that was involved in the formation of the ONC in the first place. Krumholz & Tan (2007) argued that between 6,700 and 15,000 $M_{\odot}$ were involved in the original cloud that gave rise to the ONC, and used the value of 4,500 $M_{\odot}$ quoted by Hillenbrand & Hartmann (1998) for the mass in stars. With these numbers one can compute a “present-day” efficiency between 0.3 as a minimum and 0.67 as a maximum.

We first stress that these numbers are independent on the assumed scenario of star formation. In other words, whether it occurred within one or many free-fall times, if the mass of the original cloud and of the stars are the quoted above, the present day efficiency is between 0.3 and 0.67.

It has been argued, furthermore, that while the turbulent scenario may reach such final efficiencies along many free-fall times, at a low efficiency per free-fall time, the collapsing scenario would require “extreme efficiencies” integrated over only one current, or present-day free-fall time (of the order of a few $10^{5}$ yr), of the order of $\epsilon>0.3$, in order to produce bound systems (Krumholz et al., 2019). There are several problems with this argument. The main one, as in the previous case, it assumes that the star formation occurs within one present-day free-fall time, i.e., the free-fall time estimated at the current mean gas density of the cluster. As commented out in the previous section, this time is not the relevant free-fall time involved in the hierarchical and chaotic collapse scenario, as we discussed in the previous section, but the free-fall time from the beginning of the cloud contraction, when the cloud was substantially less dense.

In addition, it should be recalled that the aforequoted 4,500$M_{\odot}$ in stars in the ONC is estimated from virial equilibrium, and, according to Hillenbrand & Hartmann (1998), $\sim$50% of that mass comes from the gas in the cloud behind the optical cluster. That reduces a factor of 2 the apparently “extreme” efficiency that the collapse model will require. But even more, recent estimations from Da Rio et al. (2012a) quote a total of 1,000 $M_{\odot}$ in stars, putting the efficiency of star formation for the ONC below 10% if the original mass of the cloud was indeed 15,000 $M_{\odot}$, as quoted by Krumholz & Tan (2007). Such efficiency could be easily achieved within one original (few Myr) free-fall timescale, without implying a “huge” star formation efficiency.

There is also the concept of “efficiency per free-fall time”, i.e., the efficiency that would have occurred after one free-fall time. Obviously, this is different from the terminal efficiency if the process of star formation has taken more than one free-fall time. Turbulence models assume that clusters are formed along $\sim$10 or more free-fall timescales, while hierarchical and chaotic collapsing models collapse within $\sim$1$-$2 free-fall timescales. Thus, the estimated efficiency per free-fall time is of the order $\sim$0.005$-$0.01 for the turbulent models, and between $0.02-0.8$ for the global collapsing model. But it should be noticed that the free-fall timescale is substantially different in each model, of the order of few $10^{5}$yr for the turbulent model, and of the order of few Myr for the hierachical and chaotic model, as we discussed in §4.2.3. The very concept of efficiency per free-fall time somehow has not much relevance if the free-fall timescale is not constant, but varies with time, as occurs in the collapse models. In these, in fact, the star formation becomes accelerated as collapse proceeds, a feature that is present in star forming regions (Hartmann et al., 2012), and that models with constant efficiency per free-fall time cannot reproduce.

Finally, it should be noticed that, part of the observational evidence that the “small” and almost independent on the density ${\rm SFR_{\rm ff}}$ is based on the assumption that the HCN traces dense gas, of the order of $6\times 10^{4}$cm^-3. However, Kauffmann et al. (2017) has shown that HCN actually traces gas with densities of the order of $n\sim$ 900 cm^-3, substantially less dense than it has been thought. In fact, if such density had been used in Fig. 5 from Krumholz & Tan (2007), it will had been difficult to argue that the ${\rm SFR_{\rm ff}}$ was small and independent on the density, as proposed by these authors.

Another wrongly posed argument against the GHCC model is that clusters have to pass through a phase in which their mass density should be larger than the mass density of the resulting stellar clusters (Krumholz et al., 2019, see arguments regarding their Fig. 14).

First of all, it should be noticed that the same problem can be posed to the turbulent model: if a volume of molecular gas forms a stellar cluster, the gas mass density of the parent cloud should be similar or larger to the stellar mass density after cluster formation and cloud dispersal, contrary to what it is shown in Krumholz et al. (2019)’s Fig. 14.

Second, this misconception probably arises because clouds are historically thought to be static. Then, if a cluster forms in the densest clump, and the final efficiency is smaller than 100%, by construction the parent clump should be more dense than the cluster it formed. But as we have argued, hierarchical and chaotic collapse is not as simple as such a picture. On the contrary, numerical simulations of collapsing clouds forming clusters actually explains with relatively simplicity the observations described by Krumholz et al. (2019). In these models, clusters form by a combination of stars forming in the central collapsing region, but also by incorporating newborn stars from the vicinity (see Figs. 1 and 2 in Kuznetsova et al., 2015). During this process, the stellar mass density can become substantially larger than the gas mass density by a combination of the inclusion of newborn stars from the vicinity, but also by gas starvation in the central region, as shown in a variety of numerical simulations (e.g., Girichidis et al., 2012; Ballesteros-Paredes et al., 2015). This process is also schematically depicted in Fig. 12.

To show that this is the case, we again look at the simulations by Kuznetsova et al. (2015, similar results are found for the other runs), where the stellar clusters are build up during the process of collapse, as gas and stars are incorporated continuously into the central cluster. In Fig. 13 we plot the mass density of gas (solid lines) and the mass density of stars (dotted lines) as a function of time, for the last 1/3 of the evolution time (recall that we evolve that simulation over one initial free-fall time). The mass densities are computed over spherical volumes of radii as indicated in the figure, and, at every time, they are centred at the position of the main cluster.

As it can be seen from Fig. 13, during the process of the cluster formation, there is a moment in which the stellar density (dotted lines) becomes larger than the gas density (solid lines). Only for large enough spheres (purple line), the gas density are always larger than the stellar densities, although it is clear that even in this case, the stellar density approaches to the gas density towards the end of the computation.

Finally, in addition to what we have discussed, it should be noticed that during the process of collapse, stellar clusters tend to develop strong density concentrations, as the entropy of gravity dominated cluster is not bound: self-gravitating systems have no maximum entropy configuration, and the more concentrated they are, the higher their entropy (Lynden-Bell & Wood, 1968, see also Binney  Tremaine, §7.3.2).

Thus, it cannot be surprising that there are no molecular clouds as dense as stellar clusters, and finding that stellar clusters do exhibit larger densities compared to the densities of molecular clouds cannot, in any way, be proof that clusters cannot be formed by a process of global collapse. In fact, rather than a problem, its a natural outcome from the GHCC model.

blackWright & Parker (2019) found evidence for an increasing velocity dispersion with mass in NGC 6530, a young ($\sim$2 Myr) stellar cluster at a distance of $\sim$1.3 kpc. They interpret this trend in velocity dispersion as a result of collapse over time to a more compact configuration where the Spitzer instability causes the massive stars to form a smaller system with a higher velocity dispersion. Wright & Parker refer to N-body simulations analyzed in Parker & Wright (2016) to argue that the observed mass-dependent velocity dispersion requires not only cool collapse but an initial highly substructured spatial distribution of stars.

blackThese observational results for NGC 6530 differ from ours for the ONC/Orion region, where we find no real evidence for a trend of velocity dispersion with stellar mass, even though we similarly argue for cluster collapse from subvirial conditions. In considering possible explanations for the difference, it is important to recognize that the velocity dispersions of both low- and high-mass stars will be strongly determined by the depth of the gravitational potential well. In the case of NGC 6530, the high mass stars tend to be much more spatially concentrated than the low-mass sample (see Figure 1 of Wright & Parker (2019)) and so consistent with higher velocity dispersions correlating with deeper gravitational potential. Perhaps more importantly, Wright et al. (2019) found evidence that on large scales the cluster is expanding, which is clearly inconsistent with a simple model of sub-virial cluster collapse. This possible expansion could be a result of expulsion of most of the initial gas in the region, as suggested by the relatively low extinction to many cluster members and the morphology of molecular gas in the region (Tothill et al., 2002), which mostly lies around the periphery of the cluster. Both of these apparent features - expansion and lack of dynamically significant molecular gas - differ qualitatively from that of our Orion/ONC sample.

blackWith regard to differences between the numerical simulations, those of Parker & Wright (2016) are pure N body simulations without gas and do not involve collapse from initially much larger scales. In addition, the simulations of Parker & Wright (2016) that show substantial growth of the velocity dispersion of the massive stars assume strongly subvirial initial conditions ($\alpha_{vir}=0.3$) and/or strong initial substructuring ($\alpha_{vir}=0.5$, fractal dimension $D=1.6$); cases of less substructuring for $\alpha_{vir}=0.5$ show little or no velocity dispersion growth. It is not obvious that our cluster simulations show similar substructuring, and this again may be affected by the presence of dynamically-important gas.