Gravity or turbulence? – IV. Collapsing cores in out-of-virial disguise

ABSTRACT

1 INTRODUCTION

1.1 The initial interpretation of the supersonic linewidths

Since the first detections of molecular gas, it was recognized that their line profiles exhibit supersonic widths (Wilson et al. 1970). Early models interpreting such profiles suggested that they were the signature of clouds in a state of large-scale collapse, since turbulence should decay over a dynamical crossing time (Goldreich & Kwan 1974; Liszt et al. 1974). This proposal was rapidly dismissed by Zuckerman & Palmer (1974), who argued that the star formation rates would be too large if clouds were in a state of free-fall. These authors proposed, instead, that the supersonic linewidths were produced by small-scale turbulence, which furthermore, might provide support to molecular clouds (MCs) against their own collapse (see e.g. Vázquez-Semadeni et al. 2000; Mac Low & Klessen 2004; Ballesteros-Paredes et al. 2007; McKee & Ostriker 2007; Klessen & Glover 2016 and references therein).

Indeed, there is a number of apparently good reasons to believe that turbulence may play an important role in the structure and dynamics of MCs. First, the Reynolds numbers in the interestellar medium are very large, indicating that the velocity field should be strongly turbulent (Elmegreen & Scalo 2004). Secondly, the velocity field in the ISM exhibits a multiscale nature that extends over several orders of magnitude in size (e.g. Larson 1981). Thirdly, clouds exhibit a fractal appearance (Falgarone, Phillips & Walker 1991), characteristic of turbulent fluids. Fourthly, stellar jets and winds, gravitational chaotic motions of orbiting gas in the Galaxy, supernova (SN) explosions, spiral arm shocks, etc., may be powering the gas of kinetic energy at multiple scales simultaneously (e.g. Norman & Ferrara 1996). In addition, the large widths of the observed spectral lines indicate that such turbulence should be supersonic, reinforcing the idea that turbulence is intrinsically related to MC structure by shaping, morphing, and fragmenting MCs (Ballesteros-Paredes, Vázquez-Semadeni & Scalo 1999a; Hopkins 2012). In the prevailing interpretation, MCs are thought to be near virial equilibrium, with perhaps a moderate excess of kinetic energy, such that some external pressure is needed to confine them for several dynamical times (Bertoldi & McKee 1992; Field, Blackman & Keto 2011; Miville-Deschênes, Murray & Lee 2016; Leroy et al. 2015). MC turbulence plays, thus, a crucial role providing global support against collapse, while promoting local collapse where motions are converging (e.g. Ballesteros-Paredes et al. 1999a; Mac Low & Klessen 2004; Ballesteros-Paredes 2006; Ballesteros-Paredes et al. 2007; Hennebelle & Falgarone 2012; Padoan et al. 2014).

1.2 Some caveats of MC turbulence

Although the turbulent nature of the velocity field in MCs is hardly questionable, several problems remain with this standard interpretation, in which turbulence can provide support against the global collapse of the clouds: First, turbulence consists of a hierarchy of coherent motions over a wide range of scales for which the largest amplitude velocity fluctuations occur at the largest scales. This implies that the dominant MC turbulent motions are far from being the small-scale turbulent motions that can provide the internal turbulent pressure to support the cloud against gravity (see e.g. Ballesteros-Paredes 2006), as confirmed numerically by Brunt, Heyer & Mac Low (2009). Thus, the dominant motions induced by turbulence at the scale of clouds will be cloud-scale motions such as compression, expansion, shearing or rotation, but not thermal-like, small-scale motions that can provide an internal pressure gradient to support it against gravity (see e.g. Ballesteros-Paredes et al. 1999a; Ballesteros-Paredes 2006; Vázquez-Semadeni 2015).

Secondly, since turbulence is a dissipative phenomenon, it requires constant driving to be maintained. The standard belief is that the turbulent energy is injected by mainly two mechanisms: one, through instabilities occurring during the assembly stage of the clouds (e.g. Vishniac 1994; Walder & Folini 2000; Koyama & Inutsuka 2002; Audit & Hennebelle 2005; Heitsch et al. 2005; Vázquez-Semadeni et al. 2006; Klessen & Hennebelle 2010). The other, by stellar feedback via outflows, winds, ionizing radiation, or SN explosions (e.g. Mac Low & Klessen 2004; Wang et al. 2010; Vázquez-Semadeni et al. 2010; Gatto et al. 2015; Padoan et al. 2016). Regarding the first point, numerical simulations of the formation of MCs from converging flows generally suggest that the turbulence level produced by various instabilities is significantly lower (⁠|${\lesssim } 20{{\rm \,per\,cent}}$|⁠) than the level typically observed in MCs (e.g. Koyama & Inutsuka 2002; Heitsch et al. 2005; Audit & Hennebelle 2005; Vázquez-Semadeni et al. 2007, 2010), even with the inclusion of external SN explosions that trigger the compressions which form the clouds (Ibáñez-Mejía et al. 2016; Seifried et al. 2018), although this is still matter of debate (e.g. Padoan et al. 2016). With respect to the stellar feedback, observations of molecular gas in nearby external galaxies does not show that the level of turbulence is correlated to the star formation rate (Caldú-Primo et al. 2013). In addition, it is not clear why the energy it injects would cause apparent near virialization of all molecular structures. In fact, numerical simulations suggest quite the opposite: clumps and moderate-mass MCs appear to be readily disrupted by photoionizing radiation or SN explosions from massive stars, while high-mass MCs are difficult to prevent from collapsing (e.g. Dale, Ercolano & Bonnell 2012; Dale, Ercolano & Bonnell 2013a; Dale et al. 2013b; Colín, Vázquez-Semadeni & Gómez 2013). Thus, neither of the two mechanisms is likely to produce the observed magnitude of the non-thermal motions in MCs.

1.3 Gravity as the driver of MC turbulence

In the last decade, numerical simulations of the formation and evolution of molecular clouds have renewed the idea that molecular clouds may be in a state of collapse. These models suggest that their formation from convergent motions in the warm neutral medium (WNM) in the Galaxy (Ballesteros-Paredes et al. 1999a; Ballesteros-Paredes, Hartmann & Vázquez-Semadeni 1999b) involves a phase transition to the cold neutral medium (CNM; Hennebelle & Pérault 1999; Koyama & Inutsuka 2000; Vázquez-Semadeni et al. 2006), and thus the clouds may be formed with masses much larger than their Jeans mass (Hartmann, Ballesteros-Paredes & Bergin 2001; Vázquez-Semadeni et al. 2007; Heitsch & Hartmann 2008; Gómez & Vázquez-Semadeni 2014), implying that they may very well be in a state of large-scale hierarchical and chaotic collapse (Vázquez-Semadeni et al. 2009; Ballesteros-Paredes et al. 2011a).

Note, however, that contrary to the suggestion by Zuckerman & Palmer (1974), such a state of global collapse does not necessarily imply that the star formation rate will be too large: the initial fragmentation occurring while the cloud is being assembled has the consequence that only very small amounts of mass are deposited in the high-density regions that are responsible for the instantaneous rate of star formation in the clouds (Ballesteros-Paredes et al. 2011b), as evidenced by the steep negative slope of the column density probability distribution functions of MCs (Kainulainen et al. 2009; Ballesteros-Paredes et al. 2011b; Kritsuk, Norman & Wagner 2011). Having large non-linear amplitudes, such dense regions complete their collapse much earlier than the whole cloud, causing a spread in the ages of the stellar products (e.g. Zamora-Avilés, Vázquez-Semadeni & Colín 2012; Hartmann, Ballesteros-Paredes & Heitsch 2012; Vázquez-Semadeni, González-Samaniego & Colín 2017). This will eventually produce massive stars that will blow strong winds and rapidly ionize the gas, effectively dispersing the nearby dense gas (e.g. Vázquez-Semadeni et al. 2010; Zamora-Avilés et al. 2012; Dale et al. 2012, 2013a, b; Colín et al. 2013; Zamora-Avilés & Vázquez-Semadeni 2014; Vázquez-Semadeni et al. 2017), shutting off the local star formation episodes by the time only |${\sim } 10{{\rm \,per\,cent}}$| of the gas has been converted to stars, thus keeping the global star formation rate and efficiency low.

In this scenario, the supersonic non-thermal motions observed in MCs, albeit still turbulent in the sense of having a chaotic component, are dominated by an infall component that occurs at multiple scales, constituting a regime of collapses within collapses (e.g. Vázquez-Semadeni et al. 2009, 2017). This implies that the non-thermal motions are dominated by a convergent component that results from the gravitational contraction, rather than being a fully random velocity field that can provide a pressure gradient capable of opposing the collapse (Vázquez-Semadeni et al. 2008; González-Samaniego et al. 2014). Nevertheless, the chaotic, moderately supersonic nature of the truly turbulent motions does produce a spectrum of density fluctuations that provides the focusing centres for the multiscale collapse (Clark & Bonnell 2005).

1.4 Evidences in favour of the hierarchical, chaotic, and gobal collapse scenario

Two important pieces of evidence supporting the scenario of global, hierarchical collapse are worth noting. On the theoretical side, as Vázquez-Semadeni et al. (2007) showed, during the early stages of the formation of a MC, the kinetic energy, driven by instabilities in the dense layer produced by the collision of diffuse gas streams, is not coupled to the gravitational energy. The two energies become correlated once the cloud becomes dominated by gravity and begins to collapse (see fig. 8 of Vázquez-Semadeni et al. 2007). In other words, the gravitational collapse is the physical agent that induces apparent virialization between the energies of the cloud.

The question of whether clouds are in a state of global, hierarchical, chaotic collapse or instead they are globally supported by turbulence against collapse for several free-fall times is crucial to the understanding of the actual dynamics of MCs, to what defines their internal structure, and to how MCs evolve and form stars. In this contribution, we present new observational and numerical evidence that massive dense cores may be in a state of collapse, even though often they may appear to have an excess of kinetic energy. In Section 2, we discuss the evolution of a contracting core in the Larson diagram and its energy budget. In Section 3, we present the observational and numerical data, while in Section 4 we show that the locus of our observational sample is similar to the locus of our collapsing cores in the simulations either in the Larson velocity dispersion-size diagram, as well as in the Larson’s ratio–column density diagram. In Section 5, we show that the apparently overvirial collapsing cores are actually not overvirial when the whole distribution of mass is taken into account on the gravitational potential. Finally, in Sections 6 and 7 we make a general discussion and provide our conclusions, respectively.

2 SCALING RELATIONS AND THE EVOLUTION OF TURBULENT CORES

2.1 Larson’s relations and virial balance

It is worth noting that the sample of Heyer et al. (2009) exhibits, at face value, systematic kinetic energy excesses with respect to the virial value. These authors argued that this might be due to the fact that the cloud masses might be systematically underestimated by a factor of ∼2. However, others have suggested instead that the excess in the Larson ratio must be indicative of the presence of an external pressure that confines the clouds (Keto & Myers 1986; Field et al. 2011; Leroy et al. 2015). We should stress, however, that pressure confinement is actually quite unlikely in the interior of molecular clouds. As extensively discussed by Vázquez-Semadeni et al. (2005), for a density fluctuation (generically, a 'clump') to attain a stable equilibrium it is necessary that the object is confined by a tenuous, warm medium that provides pressure without adding weight. This does not appear feasible in the deep interiors of MCs, where the medium is roughly isothermal, and there is no diffuse confining phase. Moreover, some of the highest column density clumps would require enormous confining pressures, ∼10⁴⁻⁹ K cm⁻³ (Field et al. 2011; Leroy et al. 2015), which appear highly unlikely to be attained by the low-density regime inside MCs.

Instead, a much simpler mechanism for clumps to exhibit a relation like equation (1) is if they are collapsing. In this case, there is no need to replenish turbulence nor to invoke pressure confinement, since, as shown by Ballesteros-Paredes et al. (2011a), the velocities produced by gravitational collapse necessarily satisfy equation (1).

However, an important question remains: although often with apparent excesses of kinetic energy (Heyer et al. 2009; Leroy et al. 2015), and sometimes with a deficiency of kinetic energy (e.g. Kauffmann, Pillai & Goldsmith 2013; Ohashi et al. 2016; Sanhueza et al. 2017; Traficante et al. 2018), why do most clouds then exhibit a correlation between this ratio and the column density? In the remainder of this contribution, we address this question.

2.2 Energy balance evolution of collapsing cores

To illustrate the expected evolution of a collapsing core in the Larson ratio–column density (⁠|${\cal L}$|–Σ) diagram, in this section we consider a spherical, uniform collapsing core that may contain both a turbulent (σ_v,3D) and a gravitationally driven (or 'infall', σ_in) components of the velocity dispersion.

It should be noted, however, that the turbulent component we consider is not necessarily assumed to provide support. Quite the contrary, as mentioned in Section 1, turbulence is known to have the largest velocities at the largest scales, and so the dominant turbulent motions in any structure must be those at the scale of the structure itself. Generally, these can be compressive, shearing, or vortical motions (Vazquez-Semadeni, Passot & Pouquet 1996; Vázquez-Semadeni et al. 2008; Federrath, Klessen & Schmidt 2008; González-Samaniego et al. 2014). Recently, Camacho et al. (2016), examining numerical simulations of the formation and collapse of MCs, have found that low column density clumps that exhibit an excess of the Larson ratio have, in roughly half of the cases, a negative average velocity divergence – i.e. a convergence. These authors interpreted this result as clumps that are being assembled by external compressive motions that are not driven by the self-gravity of the cloud, but rather constitute large-scale turbulent motions in the ambient ISM.¹ The most probable velocity gradients corresponding to convergence observed by Camacho et al. (2016) are in the range −0.5 ≲ ∇ · v ≲ 0 km s⁻¹ pc⁻¹. It is this type of non-self-gravitating motions that we have in mind when we consider turbulence in this section. We now consider how the relative contribution of these motions compare with those driven by the self-gravity of the clouds.

2.2.1 Core with no turbulence

Equation (5) implies that, as a non-turbulent core becomes unstable and begins to collapse, it follows the trajectory described by the solid line in the |${\cal L}-\Sigma$| diagram shown in the left-hand panel of Fig. 1. This can be compared to the locus of a core of the same mass and initial radius, but assuming it has the free-fall velocity at all times, given by the dotted line in this figure. This implies that cores that start at rest and have not yet attained the full free-fall speed will in general appear sub-virial by an amount that depends on the time elapsed since collapse started. Finally, note that, in reality, the actual infall speed must be even lower because, as already pointed out by Larson (1969), at the early stages of collapse of an object of mass only slightly larger than the Jeans mass, the thermal pressure is not negligible, making the collapse slower than free fall.

Figure 1.

Left-hand panel: Larson’s ratio versus column density diagram. The dark solid blue line shows the trajectory of a core of fixed mass M = 1 M_⊙ that becomes Jeans-unstable, beginning to collapse at a time t₀, at which it has an initial radius R₀ = 0.2 pc, implying an initial column density Σ₀, shown by the vertical black dash–dotted line. The black dotted line shows the locus of a core of the same mass that contracts having always the free-fall speed, as given by equation (6). The red dashed lines show the evolutionary paths of two cores with a combined turbulent+gravitational velocity, as described by equation (8), with initial velocity dispersion σ₀ = 1/2σ_ff. The three dot–one dashed light blue line corresponds to σ₀ = 3σ_ff. Of each of the last two pairs, the upper curve corresponds to η = 1/2, appropriate for a Burgers spectrum, and the lower one to η = 2, loosely representing dissipation in dense objects. Right panel: The evolution of the same set of cores (with the same labelling of the lines) in the Larson diagram, σ_v,3D versus R, represented by the dash–dotted line.

Open in new tabDownload slide

2.2.2 Core with turbulence

The right-hand panel of Fig. 1 shows the evolution of all these cases in the Larson diagram, σ_{v, 3D}−R, with the same line coding. It can be seen that, in general, the cores move transversely to the Larson linewidth–size relation, represented by the dotted curve, terminating their evolution in the upper-left part of the diagram, as first discussed by Ballesteros-Paredes et al. (2011a), and as do the cores in our numerical simulations (cf. Section 3.2 below).

Secondly, we stress that the turbulent component by no means needs to be interpreted as a supporting mechanism. As found by Camacho et al. (2016), this component corresponds, in roughly half of the low column density clouds or clumps, to large scale, external compressions that are assembling these objects, rather than supporting them.

Thirdly, the transition from a domination of turbulent to the gravitational motions is due to the fact that the two follow different scaling: while in principle the turbulent velocity dispersion depends only on scale (equation 7), regardless of column density (e.g. Padoan et al. 2016), the self-gravity-driven ones depend both on scale and column density (or mass), according to equation (1) (Heyer et al. 2009; Ballesteros-Paredes et al. 2011a; Camacho et al. 2016; Ibáñez-Mejía et al. 2016).

Fourthly, note that clumps or cores that start with a low-turbulent content may, during the initial stages of the contraction, exhibit a deficiency of the Larson ratio in relation to the equipartition (virial or free-fall) value. This may explain observations that find sub-virial cores (e.g. Kauffmann et al. 2013; Ohashi et al. 2016; Sanhueza et al. 2017; Traficante et al. 2018).

Finally, it must be remembered that in the idealized study presented here, we have considered the monolithic collapse of a single clump over two orders of magnitude in column density, from values typical of large GMCs (∼10 M_⊙ pc⁻²) to those of dense cores (∼10³ M_⊙ pc⁻²). In reality, such a range in column density is not accomplished by the collapse of a single object, but rather, through several stages of fragmentation. Thus, our result should only be taken as a first-order approximation to the effect of transitioning from an external, turbulence-dominated regime to an infall-dominated one.

3 DATA

Since we want to understand the relative importance of turbulence and/or gravity in the process of star formation, we use both observational and numerical data of massive dense cores to study the kinematics of the gas before it forms stars. To accomplish this, we need to avoid the effects of stellar feedback in both datasets. In the case of the observational data, even though our cores are located in massive star-forming regions, we choose cores without morphological evidence of perturbations (e.g. jets or free-free emission due to ultracompact H ii regions). Also, we use ammonia, which might be destroyed in case of UV stellar feedback (e.g. Palau et al. 2014), is expected to trace preferentially dense regions, rather than molecular flows, and requires large column densities to be detected. On the numerical side, stars are represented by sink particles that are allowed to accrete gas from their surrounding, but no feedback is prescribed in order to preserve the purely gravitational nature of the velocity field.

3.1 Observational data

We use recently published data of cores traced with ammonia (NH₃) and/or methylcyanide (CH₃CN), in order to map cores classified as protostellar or as quiescent and starless cores with high surface densities (N ⩾ 10²¹cm⁻²). Our sample of clumps and cores is taken from the works of Sepúlveda et al. (2011); Sánchez-Monge et al. (2013), and Hernández-Hernández et al. (2014). The idea of using this sample is to have a consistently reduced and analysed data set of different regions of intermediate to high-mass star formation, at different positions in the Galaxy, in order to avoid a particular systematic effect, as could be the case of, e.g. a wrong estimation of the distance. The first two works determine the properties from NH₃(1,1) and (2,2). The third sample includes cores studied in CH₃CN.

Special care was taken in determining the linewidths, sizes, and surface densities using the same method for all the samples. In particular, NH₃(1,1) and CH₃CN 12K-11K linewidths were measured using the same routine in the gildas program class, which allows to take into account the hyperfine structure and the blending of the lines. The NH₃ abundance was adopted to be 1 × 10⁻⁸, as an average value of previous works (Pillai et al. 2006; Foster et al. 2009; Friesen et al. 2009; Rygl et al. 2010; Chira et al. 2013), while the CH₃CN abundance has been taken as 2 × 10⁻⁹ (Hernández-Hernández et al. 2014).

It is worth noticing that, although the CH₃CN is excited at high temperatures (of the order of 100 K), the total linewidths are highly supersonic, and thus, the estimation of the non-thermal velocity dispersion is not seriously affected by the thermal broadening.

The main sources of uncertainty in our data are the abundances, the linewidths and the distances to the clouds. According to Sánchez-Monge et al. (2013), errors in the assumed abundances could lead to an uncertainty up to 60 per cent in the determination of mass, while errors in the determination of the linewidth of the order of  10–30 per cent are reasonable. Finally, the kinematic distances used in the quoted references are taken from Molinari et al. (1996), who adopt an uncertainty in the circular velocity of ±5 km s⁻¹. This is translated into an uncertainty of the order of 15 per cent in the estimation of the distance (and therefore in the estimation of the size), depending on the galactic model used (see e.g. works in Block, Freeman & Puerari 2010). However, the clumps and cores of the three samples include typical star-forming regions distributed in different locations in the Galaxy, avoiding possible biases due to uncertainties in the distance determination, or to peculiarities of a given molecular cloud. Thus, those errors are expected to affect the data in a random way, and no systematic offset is expected.

3.2 Numerical simulations

In order to numerically investigate the evolution of cores in the σ_v,3D−R ('Larson') and |${\cal L}$|–Σ diagrams, we performed five isothermal numerical simulations with Gadget-2, a Smoothed Particle Hydrodynamics (SPH) public code (Springel 2005) to represent the interior of a small (∼parsec size) molecular cloud. We also include sink particles in the code, as in Jappsen et al. (2005).

Details of the simulations can be found in Ballesteros-Paredes et al. (2015). Here, we just mention that these simulations were performed using 6 million particles, with a total mass of 1000 M_⊙ in a cubic open box of side 1 pc in all simulations but one. Three of the simulations had an initially homogeneous density field. We imposed initial velocity fluctuations using a purely rotational velocity power spectrum with random phases and amplitudes, as in Stone, Ostriker & Gammie (1998), and with a peak at wavenumber k_for = 4π/L₀, where L₀ is the linear size of the box. No forcing at later times is imposed. The initial Mach numbers of these runs were 16, 8, and 4, respectively, so we label them Runs Mach 16, Mach 8, and Mach 4, respectively.

Additionally, we performed two other runs, for which the initial velocity field is zero. The first one, labelled Run Mach 0, has initial density fluctuations with a Kolmogorov-like power spectrum proportional to k^−5/3. Finally, in the run labelled Mach 0-spherical, the initial density field is a Plummer sphere with R_c = 1 pc, a size of 5 pc and a total mass of 4500 M_⊙, again, with zero velocity field. Since the boxes are all strongly super-Jeans, all of them collapse within roughly one free-fall time. Furthermore, runs with non-zero velocity field exhibit a slight initial expansion, since we do not include any confining external medium.

The global evolution of all these runs can be generically described as follows: if the initial velocity field is not set to zero, the velocity fluctuations produce density enhancements while the turbulence is dissipated. The external parts of the cloud expand because there is no confining medium, but since the beginning, the bulk of the clouds’ mass starts to contract. This contraction, however, is noticeable only after some time, depending on the initial Mach number of the simulation: for run Mach 16, it takes about 1/2 free-fall time T_ff for the collapse to begin, ∼1/4T_ff for Mach 8 and only a small fraction of T_ff for the Mach 4 run. Once the initial imposed turbulence is dissipated, all runs proceed to collapse, lowering their sizes and increasing the velocity dispersion (although in this case, this non-thermal velocity dispersion is driven by self-gravity and corresponds to chaotic infall rather than actual turbulence that can provide support). Runs with no initial velocity field (runs Mach 0 and Mach 0-spherical), only undergo this second part of the evolution, with the initial density fluctuations driving the local centres of collapse (see e.g. Klessen & Burkert 2000).

In order to understand the evolution of these clouds in the observational diagrams (Larson ratio versus column density and velocity dispersion versus size), we have defined three regions ('cores') in the simulations for which we computed their mass, size, and velocity dispersion at every time. These are spheres located at the centre of the computational box, towards which the cloud is collapsing in a chaotic way due to the turbulent initial velocity fluctuations, and their sizes are defined so that they contain, at every time, 10, 25, and 50 per cent of the mass of the whole cloud.

The mass and the velocity dispersions are calculated straightforwardly, the former by just adding the mass of all the SPH particles sphere and the latter as the standard deviation of the particles’ velocities.

4 RESULTS

4.1 Non-existence of a σ_v,3D−R relation

As mentioned in Section 1, Ballesteros-Paredes et al. (2011a) showed that current observational datasets using high-mass MC core observations in different regions in the Galaxy do not exhibit the Larson (1981) σ_v,3D−R relation. Instead, they occupy the upper-left corner of the plane. Although such cores can be expected to have large column densities because they belong to massive star-forming regions, most of the data in the compilations discussed by Ballesteros-Paredes et al. (2011a) did not have estimations of the actual masses of the cores. Thus, the proposal that cores with larger column densities have larger velocity dispersions and smaller sizes still requires further testing using observational data that provide column density determinations that are independent of the velocity dispersion measurements.

In Fig. 2(a), we plot our observational data sample in the Larson plane σ_v,3D−R. We use σ_v,3D, defined as |$\sqrt{3} ~ \sigma _{\rm v, obs}$|⁠, assuming that the 3D velocity dispersion σ_v,3D is intrinsically isotropic, while the observed velocity dispersion σ_v,obs is only the projection of the former along the line of sight. The dataset is coloured according to the column density of each core: blue squares represent the largest column density points, ranging between 10⁵–10⁶M_⊙pc⁻²; purple triangles represent the column density range 10⁴–10⁵M_⊙; light blue, 10³–10⁴, and red diamonds, 10² and 10³M_⊙pc⁻². In addition, we highlight some cores with an asterisk, indicating those cores that are catalogued as starless by Sánchez-Monge et al. (2013).

Figure 2.

Right: Velocity dispersion versus size for the compiled observational sample. Colours and symbols denote different ranges in column density, according to the labels. Lines of constant column density if non-thermal motions are driven by self-gravity are shown as dot–dashed lines. The original fit provided by Larson (1981) is shown with the solid line. In addition, we overlap with a dark solid circle those cores that have been catalogued as quiescent and starless. Left: Same sample, in the Larson ratio |${\cal L}$| versus column density Σ space. In both frames, the open box with error bars represents the estimated uncertainties assuming an error in the estimation of the distance of 15 per cent, and considering errors of the order of 60 per cent and 30 per cent in the estimations of the mass and velocity dispersion, as estimated by (Sánchez-Monge et al. 2013).

Open in new tabDownload slide

In Fig. 2(b), we plot the ratio |${\cal L}\equiv$|σ_v,3D/R^1/2 against the column density Σ for the same sample shown in Fig. 2(a). Note that, if the clouds actually followed the two Larson relations (σ_v,3D∝R^1/2 and n∝R⁻¹) simultaneously, all data points would collapse to a single point in this plot (within the intrinsic scatter), as pointed out by Heyer et al. (2009). The red solid line in Fig. 2(b) represents the velocity dispersion of clouds in virial equilibrium, while the dashed line represents clouds in monolithic free-fall (cf. Section 2.2). We notice that, with some scatter, cores tend to follow a relation σ_v,3D/R^1/2∝Σ^1/2, although for our sample, the cores have systematically larger-than-virial values of the Larson ratio. Note that the clouds analysed by Heyer et al. (2009) also exhibit larger-than-virial values, although those authors argued that this may have been due to their clouds’ masses possibly having been underestimated by factors ∼2.

The fact that our data points are in general located above the virial-equilibrium and free-fall lines in the |${\cal L}$|–Σ plane would normally be interpreted as implying that the objects are not bound. However, this argument is somehow contradictory, since both our sample and the data of Heyer et al. (2009) contain high-column density objects, and thus they are likely to be strongly gravitationally bound. Moreover, note that, as evident in Fig. 2(b), the free-fall line lies slightly above the line for virial equilibrium; that is, the kinetic energy associated with free-fall is slightly larger than that associated with virial equilibrium. We will explore this point in Section 6.

A final comment regarding possible artefacts due to uncertainty are worth noting: as has been commented out by Sánchez-Monge et al. (2013), one could expect approximately  10–30 per cent uncertainty in the determination of the linewidth, and up to 60 per cent uncertainty in the determination of the mass of the cores, due to the line-fitting process and the chemical abundances. The size of the cores has also associated an error due to the uncertainty in the determination of the distance to the source of about 15 per cent. However, such errors should affect randomly the cores, so one should expect that the tendencies will not introduce a systematic effect that could alter the tendencies showed in this section.

4.2 Evolutionary trends of collapsing cores

In order to interpret the observational data, in particular the apparent average kinetic energy excess over equipartition of the whole sample, which is not predicted by the analytical model of Section 2.2, we now turn to the evolution of the cores in the numerical simulations discussed in Section 3.2.

In Fig. 3, we show the evolution of the simulated cores in the Larson plane. As in Fig. 2(a), the dotted line represents Larson’s original relation, with a slope of 0.38, while the dot–dashed lines represent lines of constant column density. For reference, the yellow area denotes the region occupied by the observed cores shown in Fig. 2(a). In each panel, the evolution of the numerical cores is indicated with dashed lines, with red, dark purple, and cyan representing the cores containing 10, 20, and 50 per cent of the total mass in the box, respectively. The solid lines are discussed in the next section. The length of these trajectories corresponds to one free-fall time at the initial mean density for runs Mach 16, Mach 8, and Mach 4, and to 0.8 free-fall times for the zero-velocity runs (Mach 0 and Mach 0-spherical), since these simulations contain very dense regions and then the time-step becomes too small, slowing the simulation down towards the end of the evolution. The triangles at the extreme of these lines indicate the place in which each core terminates its evolution at the end of the simulation. In all runs, the cores with the higher masses have larger radii, so that the rightmost curves represent the evolution of the core with 50 per cent of the mass, while the leftmost ones correspond to those with only 10 per cent of the mass.

Figure 3.

Dashed lines: Evolution of the cores in the simulations in the Larson plane. Panels are organized as follows: (a) run with initial Mach number 16, (b) Mach 8, (c) Mach 4, (d) Mach 0 with Kolmogorov density fluctuations, and (e) Mach 0, spherical with a Plummer profile. Each curve represents a region that contains 50 per cent (rightmost curve), 25 per cent, and 10 per cent (leftmost curve) of the mass at each time. Dashed lines are lines where energy equipartition is achieved at different column densities (see equation 6). Solid lines: Evolution of the cores in the simulations, but considering the correction factor defined in Section 5. Triangles and squares denote the last time in the evolution.

Open in new tabDownload slide

The evolution of the cores in the Larson diagram can be described as follows: Runs Mach 16 and Mach 8 initially exhibit a decrease in the velocity dispersion. The former does so at nearly constant size, while the latter does so while slightly lowering its size. This decrease corresponds to the initial period during which the excess turbulent energy dissipates. But, as gravity takes over, the sizes start decreasing more rapidly and the velocity dispersion starts increasing again, causing the cores to move towards the upper left in the Larson diagram. In the case of runs with small initial turbulence (Mach 4) or without it (Mach 0 and Mach 0-spherical), the cores only undergo the second stage, lowering their size while increasing their velocity dispersion.

_v,3DIn Fig. 4, we show the Larson ratio |${\cal L}=$|σ/R^1/2 versus the mass column density Σ of the cores in our simulations. The lower and upper red dotted lines with slope 1/2 represent the loci of a sphere with uniform density and (a) |E_g| = 2E_k (virial equilibrium) and (b) |E_g| = E_k (energy equipartition, expected for collapse), respectively. The yellow shaded area again denotes the region occupied by our observational data (see Fig. 2b).

Figure 4.

Dashed lines: Evolution of the cores in the simulations in the Larson ratio versus column density plane. Panels and notation, as in Fig. 3. The two lines in this figure denote the locus of the virial and energy equipartition of a sphere with constant density.

Open in new tabDownload slide

The evolution of our cores in the Larson ratio–column density Σ diagram is represented, again, by the dashed lines: red, dark purple, and cyan representing the cores containing 10, 20, and 50 per cent of the total mass in the box, respectively. As in the previous figure, the solid lines are discussed in the next section. As can be expected, runs with an excess of initial kinetic energy (Mach 16, Mach 8) lower their |${\cal L}$| ratio during the initial dissipation period, but as soon as gravity takes over, they increase it again, while their column density also increases. Instead, runs with low initial turbulence (Mach 4) or without it (Mach 0 and Mach 0-spherical) increase their Larson ratio |${\cal L}$| while they also increase their column density, in agreement with the prediction of the analytical model.

What seems to be surprising in Fig. 4 is that, even though the cores we show are all collapsing, they terminate their evolution exhibiting Larson ratios in excess of the equipartition and virial values corresponding to their column density. Thus, if only the information of this plot were provided, they would be naively interpreted either as expanding cores, or as cores that require an external confining pressure in order to be in virial equilibrium and to avoid their expansion.

Another point to notice is that the final column densities of the cores that would be inferred from their final position in the Larson diagram (comparing with the constant-column density lines in Fig. 3) are larger than their actual final column densities (as denoted by the abscissa of the triangles in Fig. 4). For instance, the Larson diagram would imply a column density of the order of ∼10⁴ M_⊙pc⁻² for the cores from run Mach 16 (upper-left panel in Fig. 3), although they actually end their evolution with a column density of Σ ∼ 10³M_⊙pc⁻² (see the upper-left panel in Fig. 4). In the following section, we will show that these results can be easily explained when the true gravitational energy of the cores, i.e. the gravitational energies considering the complex internal density structure.

5 COLLAPSING CORES IN SUPERVIRIAL DISGUISE: NEED FOR A CORRECTION FACTOR

It is clear from the previous section that observed cores do not follow Larson’s velocity dispersion–size scaling relation, and may exhibit apparent excesses of kinetic energy when compared to the gravitational energy of a spherical cloud with constant density of the same mass and size, as is implicitly done by the standard virial equilibrium and equipartition lines drawn in Figs 2(b) and 4.

In the standard picture, clouds/cores with such excesses are interpreted either as expanding (e.g. Dobbs, Burkert & Pringle 2011), or as confined by an external medium (e.g. Keto & Myers 1986; Bertoldi & McKee 1992; Field et al. 2011). A frequent explanation of the origin of this apparent kinetic energy excess in those cores is that they belong to massive star-forming regions, and thus they may be subject to kinetic energy injection from stellar feedback. However, at least in the case of the cores in our observational sample, the observed molecules, ammonia (NH₃) and methylcyanide (CH₃CN), are thought to not be seriously affected by stellar feedback: Ammonia is destroyed by radiation and winds from stars, while CH₃CN requires large column densities to be detected, and thus it is unlikely that it comes from the outflows. In addition, even those cores that are known to be starless and quiescent in our sample (denoted by overlaid asterisks in Fig. 2b) also exhibit an excess of kinetic energy. How this apparent excess should be interpreted?

We have also calculated the correction factor for the three cores in our observational sample that have reported density profiles. Details on the method for the calculation of the density profiles can be found in Palau et al. (2014). Here, we just mention that the radiative transfer equation is solved assuming that the density and temperature profiles have power-law dependencies in the range 0 ⩽ r ⩽ 10⁴au, and then comparing the observed intensity profile in the sub- and millimetric regimes, as well as the resulting spectral energy distribution, to the observational data. From the 19 sources reported in Palau et al. (2014), 5 cores are present in our sample, but 2 of them have power-law indices ⩾2, implying that their gravitational energy diverges.

In Fig. 5, we show the velocity dispersion–size relation and the Larson’s ratio versus column density of these cores. As it could be expected, the Larson’s ratio after correction is smaller, and thus, cores that appear to be overvirial now become sub- or near-virial. This shows again the need for better estimations of the gravitational content when trying to infer the actual gravitational state of cores, and the role of turbulence in supporting clouds.

Figure 5.

Velocity dispersion versus size (left), and Larson’s ratio versus column density (right) for three of the cores in our sample for which Palau et al. (2014) have fitted their density structure, which corresponds to power laws with indices flatter than −2. Blue points are the data as plotted in Fig. 2, while green points represent the data after correcting the y axis by the Γ factor. As it could be expected, the green points exhibit smaller values of the velocity dispersion (left-hand panel) and Larson’s ratio (right-hand panel), showing that not necessarily the cores are overvirialized when they appear at first glance above the ‘virial' or ‘collapse' lines. This reinforces the need for better determinations of the density structure in order to understand the actual dynamical state of molecular cloud cores.

Open in new tabDownload slide

6 DISCUSSION

Traditionally, the apparent near-virialization exhibited by clouds and their substructures (clumps and dense cores) has been interpreted as a manifestation that all the structures in this hierarchy are supported against their self-gravity by strongly supersonic turbulence, and that, when this turbulence dissipates at the smallest scales, the cores can then proceed to collapse (see e.g. the reviews by Larson 1981; Vázquez-Semadeni et al. 2000; Mac Low & Klessen 2004; Ballesteros-Paredes et al. 2007; McKee & Ostriker 2007; Bergin & Tafalla 2007; Donkov, Veltchev & Klessen 2011, 2012; Hennebelle & Falgarone 2012; Dobbs et al. 2014; Veltchev, Donkov & Klessen 2016). However, it is difficult to understand how the many energy-injection mechanisms could adjust themselves to provide just the right amount of energy to the turbulence to keep the structures nearly virialized at all scales. Some authors have proposed, on the basis of idealized models for the star formation rate (e.g. Krumholz, Matzner & McKee 2006; Goldbaum et al. 2011), that the feedback from stellar sources internal to the clouds can self-regulate to attain near virialization. However, evolutionary models (Zamora-Avilés et al. 2012; Zamora-Avilés & Vázquez-Semadeni 2014), as well as numerical simulations (e.g. Vázquez-Semadeni et al. 2010; Dale et al. 2012, 2013a,b; Colín et al. 2013) do not show a trend towards local virialization, but rather, towards local destruction of the star-forming sites.

On the other hand, it is also sometimes suggested that feedback from sources external to the clouds can drive the turbulence at all scales within the clouds (e.g. Padoan et al. 2016). The underlying assumption here is that this is a natural consequence of the turbulent cascade, as the energy spectrum of Burgers-like strongly supersonic turbulence should approach the form E(k)∝k⁻², where k is the wavenumber. This spectrum implies a velocity dispersion–size scaling relation of the form σ∝R^1/2, similar to Larson’s scaling, and thus strongly supersonic turbulence has been proposed as the origin of the Larson scaling. However, as shown in several observational and numerical studies (Heyer et al. 2009; Ballesteros-Paredes et al. 2011a; Leroy et al. 2015; Camacho et al. 2016; Traficante et al. 2018 this work; although see Padoan et al. 2016), the Larson scaling is not universal, and instead the Larson ratio |${\cal L}$|⁠, which should be constant for strongly supersonic turbulence, is actually dependent on the column density as Σ^1/2, suggesting that the process does not rely only on the turbulent cascade.

More specifically, within the scatter, cores tend to be organized by column density along lines with slope 1/2 in the Larson σv,3D–R diagram, with high column density cores being located on the upper-left part of the diagram, and low column density cores appearing at lower velocity dispersions at a given value of R. Moreover, regardless of the observational certainties, not all cores with the same column density will necessarily be located exactly at the same position in the |${\cal L}$|–Σ diagram, since their exact position on the diagrams depends on their evolutionary stage, as the results from Sections 2.2 and 4.2 show.

It is worth nothing that most of the low column density cores showed by Camacho et al. (2016) in numerical simulations exhibit large Larson’s ratios. If the Γ correction were estimated for the low column density cores of Camacho et al. (2016), it is likely that a substantial population of them will appear as undervirial, or virial. In fact, half of those cores do have negative mean values of the divergency, suggesting that indeed, those are wrongly estimated as supervirial.

The observational sample considered here also shows a correlation between the Larson ratio |${\cal L}$| and the column density in general agreement with equation (1), although on average it exhibits somewhat supervirial values when the gravitational energy of the sample cores is computed using equation (11) (e.g. Heyer et al. 2009; Leroy et al. 2015), which assumes that the cloud is a uniform-density sphere. Although we cannot discard stellar feedback in the observations, the facts that (a) our sample was observed with tracers that should not be strongly affected by feedback, (b) the subsample of quiescent starless cores exhibits the same behaviour, and (c) simulations of collapsing cores also exhibit overvirialization, lead us to conclude that such apparent overvirialization may be spurious. Indeed, we find that, when the kinetic energy is corrected by a factor equal to the ratio of the actual gravitational energy to that given by equation (11), then the apparent excess of kinetic energy essentially disappears.

On the other hand, surveys of massive star-forming cores sometimes exhibit subvirial values of the Larson ratio, a result which has often been interpreted as the cores being supported by some form of energy other than the turbulent pressure, such as the magnetic pressure (e.g. Kauffmann et al. 2013; Ohashi et al. 2016; Sanhueza et al. 2017), or else as the cores being already in full-blown collapse, i.e. at the full free-fall speed, rather than the lower infall speed given by equation (4). However, because the infall speeds are also non-thermal motions, full-blown collapse should correspond to a virial ratio (cf. equation 9) α ∼ 2, not to α < 1. Although magnetic support for these cores is indeed a plausible explanation, we have shown that, when a core begins its own local collapse, it may start with subvirial velocities if its inertial turbulent energy is sufficiently low, approaching the virial value from below. This constitutes an alternative plausible explanation for the observations of subvirial cores.

The notion that the non-thermal motions observed in MCs are produced by the collapse itself is contrary to the widespread notion that turbulence is the main physical process controlling the internal dynamics of MCs, providing support to MCs and causing their fragmentation. Although the turbulent density fluctuations in the diffuse medium may very well play a crucial role in the formation of the seeds of what eventually will grow as cores via instabilities, the fluctuations produced self-consistently by this mechanism are generally not strong enough to become locally Jeans unstable, until global collapse has significantly reduced the mean Jeans mass in the cloud (e.g. Koyama & Inutsuka 2002; Clark & Bonnell 2005; Heitsch et al. 2005; Vázquez-Semadeni et al. 2007; Heitsch & Hartmann 2008). Moreover, observational data show that fragmentation levels within dense cores do not correlate with the intensity of the observed non-thermal motions (Palau et al. 2015), again suggesting that they actually do not consist of random turbulent fluctuations. Nevertheless, the fluctuations produced by the growth of seeds triggered by the global collapse do contain a significant fraction of chaotic motions due to the turbulent background, and so these fluctuations are far from being ordered and monolithic. Indeed, Heitsch, Ballesteros-Paredes & Hartmann (2009) showed that clouds in a state of hierarchical gravitational collapse do exhibit line profiles similar to observed MCs in the structure of the CO lines, the supersonic widths, and their core-to-core velocity dispersion. However, the net average component of the velocity field continues to be dominated by global contraction, as shown by studies of the dense regions in simulations of driven, isothermal turbulence, which indicate that the overdensities tend to have a negative net velocity divergence (i.e. a convergence) (e.g. Vázquez-Semadeni et al. 2008; González-Samaniego et al. 2014; Camacho et al. 2016). Thus, the density fluctuations are in general contracting, rather than being completely random with zero or positive net divergence, as would be necessary for the bulk motions to exert a 'turbulent pressure' capable of opposing the self-gravity of the overdensities. If the non-thermal motions in the clumps and cores do not exert a turbulent pressure capable of providing support against the self-gravity of the structures, then MC models based on the competition between gravity and turbulent support may need to be revisited (e.g. McKee & Tan 2003; Krumholz & McKee 2005; Hennebelle & Chabrier 2008, 2011; Hopkins 2012).

Finally, it is worth mentioning that the behaviour predicted by the models in Section 2, and reproduced by the evolutionary trend of our cores in the simulations suggests that this is a physical result. We do not expect that local variations in the abundances of the tracer molecules may modify substantially the trend, as discussed in Section 4.1

7 CONCLUSIONS

In this paper, we have used observational data and numerical simulations to show that the scaling relation between velocity dispersion and size does not hold when column densities spanning a large dynamic range are considered. In this case, cores with large column densities tend to be located in the upper-left corner of the Larson velocity dispersion–size diagram. Using numerical simulations, we showed that, as cores collapse, their sizes become smaller and their column densities larger, and thus, their evolution tends to follow lines oblique to the σ_v,3D∝R^1/2 relation.

Additionally, we showed analytically that, as cores evolve from when they first detach from the global flow and begin their local collapse, they may exhibit subvirial velocities if their internal turbulent component is initially low enough. This is because, when the core first starts to collapse locally, its gravitationally driven speed starts out from zero, and takes a finite amount of time to reach the full free-fall speed. This behaviour was observed also in the numerical simulations that start with little or no turbulent energy, and it may explain recent observations of apparently sub-virial cores (e.g. Kauffmann et al. 2013; Ohashi et al. 2016; Sanhueza et al. 2017; Traficante et al. 2018).

We also showed that, although the observed cores were selected to avoid stellar feedback, they appear to be supervirial. However, this feature is likely due to the fact that, rather than comparing the kinetic energy to the actual gravitational content of the core, W, in practice it is customary to use the gravitational energy of a sphere with constant density as a proxy for W. This proxy actually provides only a lower limit (in absolute value) to the actual value of W. We show that, indeed, this is the case in our simulations: collapsing cores do exhibit overvirial velocities at the end of their evolution, and thus end up in the virial/equipartition zone. Thus, the virial parameter computed using the uniform sphere approximation should be taken with caution when interpreting the observational data. This result is in agreement with the recently observational study by Traficante et al. (2018), who find that the value of the α_vir parameter does not correlate with the presence or not of infall signatures in cores.

This work then provides support to the scenario that non-thermal motions in MCs have largely a gravitational origin and are dominated by infall motions. Because of the rapid dissipation of turbulence, the conversion of the infall motions into random, truly turbulent motions that can oppose the collapse that produces them in the first place does not appear feasible. In a future contribution, we plan to further explore this possibility.

ACKNOWLEDGEMENTS

This work was supported by UNAM-PAPIIT grant number IN110816 to JBP and CONACYT grant 255295 to EVS. In addition, JBP acknowledges the hospitality of the Institute for Theoretical Astrophysics of the University of Heidelberg, as well as UNAM’s DGAPA-PASPA Sabbatical program. He also is indebted to the Alexander von Humboldt Stiftung for its valuable support. RSK acknowledges support from the Deutsche Forschungsgemeinschaft in the Collaborative Research Center SFB 881 ‘The Milky Way System' (subprojects B1, B2, and B8) and in the Priority Program SPP 1573 ‘Physics of the Interstellar Medium’ (grant numbers KL 1358/18.1, KL 1358/19.2). RSK furthermore thanks the European Research Council for funding in the ERC Advanced Grant STARLIGHT (project number 339177). Numerical simulations were performed in the supercomputer Miztli, at DGTIC-UNAM. We have made extensive use of the NASA-ADS database.

Footnotes

REFERENCES