Global Hierarchical Collapse In Molecular Clouds. Towards a Comprehensive Scenario

TL;DR

This paper presents a comprehensive model of Global Hierarchical Collapse in molecular clouds, explaining star formation, filament dynamics, and mass distribution through a unified collapse scenario.

Contribution

It introduces a unified GHC framework that describes multi-scale collapse, filament formation, and star formation processes in molecular clouds, integrating various observational features.

Findings

Filaments accrete onto hubs, forming a collapse flow from cloud to core scales.

Molecular and dense mass fractions increase over time in clouds.

Prestellar cores are contracting and may appear unbound early on.

Abstract

We present a unified description of the scenario of Global Hierarchical Collapse and fragmentation (GHC) in molecular clouds (MCs), owing to the continuous decrease of the average Jeans mass in the contracting cloud. GHC constitutes a regime of collapses within collapses, in which small-scale collapses begin at later times, but occur on shorter timescales than large-scale ones. The difference in timescales allows for most of the clouds' mass to be dispersed by feedback from the first massive stars, maintaining the global star formation rate low. All scales accrete from their parent structures. The main features of GHC are: star-forming MCs are in an essentially pressureless regime, which produces filaments that accrete onto clumps and cores ("hubs"). The filaments constitute the collapse flow from cloud to hub scales and may approach a quasi-stationary state; the molecular and dense mass fractions of the clouds increase over time; the first (low-mass) stars appear several Myr after global contraction began; more massive stars appear after a few Myr in massive hubs resulting from the collapse of larger scales; the minimum fragment mass may extend well into the brown-dwarf regime; Bondi-Hoyle-Lyttleton accretion occurs at the protostellar and core scales, accounting for a near-Salpeter IMF; the extreme anisotropy of the filamentary network explains the difficulty in detecting large-scale infall signatures; the balance between inertial and gravitationally-driven motions in clumps evolves during the contraction; prestellar cores adopt Bonnor-Ebert-like profiles, but are contracting ever since early times when they may appear to be unbound and to require pressure confinement; stellar clusters develop radial age and mass segregation gradients. Finally, we discuss the incompatibility between supersonic turbulence and the observed scalings in the molecular hierarchy.

Tables

Table 1. Table 1: Comparison of the magnetic support, GT, and GHC scenarios

MC feature	Magnetic supporta	Gravoturbulent (GT)b	GHC
Time dependence	Quasi-stationary	Quasi-stationary	Evolutionary at cloud scale
Cloud formation	…	…	Large-scale compressions due to turbulence or large-scale gravity (stellar spiral arms, magneto-Jeans, etc.)
Cloud-scale supporting pressure	Magnetic	Turbulent	None
Scale of gravitational contraction	Whole cloud (for magnetically supercritical clouds). Dense core (for subcritical clouds)	Dense core	Whole cloud (possibly including cold atomic envelope)
Filament formation and structure	…	Shocks. No implication of longitudinal flow	Anisotropic gravitational contraction with longitudinal flow (river-like structures running down the gravitational potential)
Core formation and collapse mechanism	Slow gravitational contraction mediated by ambipolar diffusion	Shock compression and local reduction of the Jeans mass	Sequential destabilization of successively smaller mass scales as global Jeans mass decreases
Distinction between low- and high-mass star-forming regions	Subcritical vs. supercritical clouds	…	Evolution from low- to high-mass regions
SFR	Stationary. Low (high) in subcritical (supercritical) clouds	Stationary	First increasing due to cloud contraction, then decreasing due to feedback, over timespan $\sim 10$ Myr
IMF	…	Determined by turbulent core mass function	Bondi-Hoyle-Lyttleton-type accretion can account for $-2$ slope of power-law part of stellar IMF and core mass function. Massive stars form during peak SFR
Reason for dismissal	Most clouds are supercriticalc	(Proposed) Inconsistencies between turbulence and observed cloud properties; effect of turbulence is to form/destroy clouds, not support. Valid on average over large spatial or temporal scales	TBD

Equations

$$\frac{{\sigma }_v}{R^{1/2}}\approx \sqrt{G\Sigma },$$

$$M_{\mathrm{J}}=\frac{{\pi }^{5/2}}{6}\frac{c_{\mathrm{s}}^3}{G^{3/2}{\rho }^{1/2}}.$$

$$L_{\mathrm{J}}={\left(\frac{\pi c_{\mathrm{s}}^2}{G\rho }\right)}^{1/2}$$

$$v_{\mathrm{ff}}=\sqrt{\frac{2GM}{R}},$$

$$v_{\mathrm{g}}=\sqrt{2GM\left(\frac{1}{R}-\frac{1}{R_0}\right)}.$$

$${\mathcal{L}}_{\mathrm{g}}=\sqrt{2\pi G\Sigma \left[1-{\left(\frac{{\Sigma }_0}{\Sigma }\right)}^{1/2}\right]},$$

$${\sigma }_{\mathrm{turb}}=v_0{\left(\frac{R}{R_0}\right)}^{\eta },$$

$${\mathcal{L}}_{\mathrm{tot}}=\frac{{\sigma }_{\mathrm{tot}}}{R^{1/2}}=\sqrt{2\pi G\Sigma \left[1-{\left(\frac{{\Sigma }_0}{\Sigma }\right)}^{1/2}\right]+\frac{{\sigma }_{\mathrm{turb}}^2}{R}},$$

$${\alpha }_{\mathrm{vir}}$$

$$R(\tau )=R_0(1-{\tau }^2{)}^{a/3},$$

$$\rho (\tau )={\rho }_0(1-{\tau }^2{)}^{-a},$$

$$M_{\mathrm{J}}(\tau )=M_{\mathrm{J}}({\rho }_0)(1-{\tau }^2{)}^{a/2}.$$

$${\rho }_{\mathrm{rms}}(\tau )\sim {\mathcal{M}}_{\mathrm{s}}^2\rho (\tau ),$$

$$M_{\mathrm{J},\mathrm{rms}}(\tau )\approx \frac{M_{\mathrm{J}}({\rho }_0)}{{\mathcal{M}}_{\mathrm{s}}}{\left(1-{\tau }^2\right)}^{a/2}.$$

$${\tau }_{\mathrm{M},\mathrm{u}}\equiv \frac{t_{\mathrm{M},\mathrm{u}}}{t_{\mathrm{ff},0}}\approx {\left[1-{\mu }^{2/a}\right]}^{1/2},$$

$$M_{\mathrm{J},n}=\frac{M_{\mathrm{J},n-1}}{n_{\mathrm{MJ}}}.$$

$$M_{\mathrm{J},\mathrm{rms},1}=\frac{M_{\mathrm{J},\mathrm{rms},0}}{n_{\mathrm{MJ}}}=\frac{1}{n_{\mathrm{MJ}}}\frac{M_{\mathrm{J},0}}{{\mathcal{M}}_{\mathrm{s}}},$$

$$M_{\mathrm{J},\mathrm{rms},2}$$

$$M_{\mathrm{J},\mathrm{rms},n}=\frac{M_{\mathrm{J},\mathrm{rms},n-1}}{n_{\mathrm{MJ}}}=\frac{1}{n_{\mathrm{MJ}}^n}\frac{M_{\mathrm{J},0}}{{\mathcal{M}}_{\mathrm{s}}}.$$

$$M_{\mathrm{J},\mathrm{rms}}({\tau }_n)\equiv M_{\mathrm{J},\mathrm{rms},n}=\frac{M_{\mathrm{J},0}}{{\mathcal{M}}_{\mathrm{s}}}{\left(1-{\tau }_n^2\right)}^{a/2}.$$

$${\tau }_n=\sqrt{1-{\left(\frac{1}{n_{\mathrm{MJ}}}\right)}^{2n/a}},$$

$$M(n\ge n_{\mathrm{SF}})\ge M_{\mathrm{J}}(n_{\mathrm{SF}}),$$

$$P(s)=\frac{1}{\sqrt{2\pi {\sigma }_s^2}}\exp \left[-\frac{(s-⟨s⟩{)}^2}{2{\sigma }_s^2}\right],$$

$$M(n\ge n_{\mathrm{SF}})$$

$${\omega }^2=c_{\mathrm{s}}^2{k}^2-4\pi G{\rho }_0,$$

$$k_{\mathrm{J}}\equiv {\left(\frac{4\pi G{\rho }_0}{c_{\mathrm{s}}^2}\right)}^{1/2}$$

$$\mbox{Inner core:}\left\{\begin{array}[]{ll}\rho(r)={\mbox{c}st.}\\ v(r)\propto-r.\end{array}\right.$$

$$\mbox{Envelope:}\left\{\begin{array}[]{ll}\rho(r)\propto r^{-2},\\ v(r)=\mbox{cst.}<0.\\ \end{array}\right.$$

$$\mbox{Inner core:}\left\{\begin{array}[]{ll}\rho(r)\propto r^{-3/2},\\ v(r)\propto-r^{-1/2}.\end{array}\right.$$

$$\mbox{Envelope:}\left\{\begin{array}[]{ll}\rho(r)\propto r^{-2},\\ v(r)=\mbox{cst.}<0.\\ \end{array}\right.$$

Full text

††pagerange: Global Hierarchical Collapse In Molecular Clouds. Towards a Comprehensive Scenario–Global Hierarchical Collapse In Molecular Clouds. Towards a Comprehensive Scenario

1 Introduction

Molecular clouds (MCs) have long been known to exhibit supersonic linewidths (Wilson et al., 1970) and to possess masses much larger than their thermal Jeans mass (e.g., Blitz, 1993). One of the first interpretations of the supersonic linewidths observed in MCs (Wilson et al., 1970) was that the clouds are in global gravitational collapse (Goldreich & Kwan, 1974), but this suggestion was quickly rejected by Zuckerman & Palmer (1974) and Zuckerman & Evans (1974) on the basis of two main arguments: One, that, if the clouds were in global collapse, then the star formation rate (SFR) would be much larger than the actual observed value. Second, radial, cloud-scale motions should cause absorption lines produced by the gas surrounding star-forming regions to be systematically redshifted from emission lines produced at the region itself, since the surrounding gas should be infalling into the star-forming region. Since this systematic shift had not been observed (but see Barnes et al., 2018, for a possible recent detection), Zuckerman & Palmer (1974) ruled out global collapse, as well as any cloud-scale systematic, radial motions. As a plausible alternative, Zuckerman & Evans (1974) suggested that the observed supersonic linewidths in MCs correspond instead to small-scale turbulence, so that it can provide a roughly isotropic pressure gradient. This view has prevailed until today.

In the present-day scenario (see, e.g., the reviews by Mac Low & Klessen, 2004; Ballesteros-Paredes et al., 2007; McKee & Ostriker, 2007; Bergin & Tafalla, 2007; Hennebelle & Falgarone, 2012), the clumps and cores within MCs are thought to be produced by the supersonic compressions caused by the turbulence in the cloud, and the turbulence is expected to cascade from the large scales where it is injected all the way down to the smallest (dense core, $\la 0.1$ pc) scales, so that the clumps at intermediate scales are themselves internally turbulent. In turn, the clumps are believed to be supported against their self-gravity by their internal turbulence, and to collapse once the latter has been dissipated, so that the turbulent support is lost (e.g., Goodman et al., 1998; Bergin & Tafalla, 2007), or when the compression is strong enough that the local Jeans mass becomes smaller than the clump’s own mass (e.g., Sasao, 1973; Padoan, 1995; Vázquez-Semadeni et al., 2003a; Gómez et al., 2007; Gong & Ostriker, 2009).

However, over the last decade, a number of recent observational studies have suggested that star-forming MCs may be in a state of global gravitational contraction after all. On the basis of comparisons between millimeter interferometric observations and numerical simulations, Peretto, Hennebelle & André (2007) suggested that the elongated clump NGC 2264-C may be in the process of collapsing and fragmenting along its long axis. Several molecular-line kinematic studies of clouds and their internal filaments suggest that these systems are undergoing global multi-scale and multi-epoch collapse (e.g., Galván-Madrid et al., 2009; Peretto et al., 2013; Beuther et al., 2015; Liu et al., 2015, 2016b; Friesen et al., 2016; Jin et al., 2016; Hacar et al., 2017; Csengeri et al., 2017; Yuan et al., 2018; Jackson et al., 2019; Barnes et al., 2019), with the filaments in particular exhibiting longitudinal flow (e.g., Sugitani et al., 2011; Kirk et al., 2013; Fernández-López et al., 2014; Motte et al., 2014; Peretto et al., 2014; Tackenberg et al., 2014; Jiménez-Serra et al., 2014; Hajigholi et al., 2016; Wyrowski et al., 2016; Juárez et al., 2017; Rayner et al., 2017; Lu et al., 2018; Baug et al., 2018; Gong et al., 2018; Ryabukhina et al., 2018; Dutta et al., 2018; Chen et al., 2019) that feeds central “hubs”, the sites where the filaments converge (Myers, 2009).

On the numerical side, self-consistent simulations of giant molecular cloud (GMCs) formation by converging motions in the warm diffuse atomic medium, as suggested by Ballesteros-Paredes et al. (1999b), Hennebelle & Pérault (1999), and Hartmann, Ballesteros-Paredes, & Bergin (2001), have shown that

The clouds are dynamical entities that accrete from their diffuse environment and therefore grow in mass (e.g., Ballesteros-Paredes et al., 1999a; Audit & Hennebelle, 2005; Banerjee et al., 2009). 2. 2.

The turbulence generated by inherent instabilities in the compressed material (e.g., Koyama & Inutsuka, 2002; Audit & Hennebelle, 2005; Heitsch et al., 2005, 2006; Vázquez-Semadeni et al., 2006; Wareing et al., 2019) is only moderately supersonic, with typical sonic Mach numbers ${\cal M}_{\rm s}\sim 3$. This value corresponds to that observed at the scale of individual “molecular clumps” ($\la 1$ pc), but falls significantly short of the typical values for large clouds and GMCs (see, e.g., Heyer & Brunt, 2004). 3. 3.

Because of this relatively low Mach number, the turbulence is not enough to support a GMC, and so the cloud begins to undergo global collapse shortly after its mass becomes larger than its thermal Jeans mass (e.g., Vázquez-Semadeni et al., 2007, 2010; Heitsch et al., 2008b; Carroll-Nellenback et al., 2014). 4. 4.

The collapse process, however, is extremely nonlinear, being initially very slow, and violently accelerating towards its final stages (Girichidis et al., 2014; Zamora-Avilés & Vázquez-Semadeni, 2014). 5. 5.

Nevertheless, no star formation occurs during the first several megayears of the global collapse, both in the non-magnetic (e,g, Vázquez-Semadeni et al., 2007; Heitsch & Hartmann, 2008; Carroll-Nellenback et al., 2014) and the magnetized case (e.g. Vázquez-Semadeni et al., 2011; Fogerty et al., 2016; Zamora-Avilés et al., 2018), due to the slowliness of the initial stages of the collapse. 6. 6.

As shown by Camacho et al. (2016), Ibáñez-Mejía et al. (2016), and Ballesteros-Paredes et al. (2018), the clouds and clumps defined in the simulations follow the scaling relation observed for molecular clouds,

[TABLE]

(e.g., Field et al., 2011; Heyer et al., 2009; Leroy et al., 2015), where $\sigma_{v}$ is the cloud velocity dispersion, $R$ is its approximate radius, and $\Sigma$ is its column density. On the other hand, Ibáñez-Mejía et al. (2016) have shown that this scaling is not reproduced in simulations of the ISM with supernova-driven turbulence when the self-gravity is not included. Instead, purely supernova-driven turbulence in the ISM generates a scaling of the form $\sigma_{v}\propto R^{1/2}$ (i.e., $\sigma_{v}R^{-1/2}=$ cst., as corresponds to strongly supersonic turbulence; e.g., Padoan et al., 2016), which corresponds mainly to the scaling observed in GMCs in the outer parts of galaxies, but is not representative of the majority of star-forming MCs (e.g., Heyer et al., 2009; Ballesteros-Paredes et al., 2011a; Leroy et al., 2015; Traficante et al., 2018a).

Item (vi) above suggests that self-gravity may be the main driver of the motions at all scales in a cloud, and not necessarily implying virialization as it is frequently interpreted, but rather infall, as will be discussed in Sec. 6.4.

In addition, Burkert & Hartmann (2004) investigated the collapse of thin, finite sheets of cold gas, formed by converging flows in the warm atomic gas. On the basis of those results, Hartmann & Burkert (2007) suggested that the morphology of the Orion A cloud, and the location of the Orion Nebula Cluster within it, may be explained simply in terms of a global gravitational contraction of the cloud. Recently, evidence for this cloud-scale contraction in Orion has been presented by Hacar et al. (2017).

In the light of this mounting evidence of global MC gravitational contraction, in Vázquez-Semadeni et al. (2009, see also Vázquez-Semadeni et al. 2017), we suggested that MCs in general might be undergoing a process of global, hierarchical collapse (GHC), in which the clouds develop multi-scale gravitational contraction, similar to the hierarchical fragmentation scenario proposed by Hoyle (1953), but seeded with nonlinear density fluctuations produced by the turbulence. This mechanism constitutes a mass cascade, somewhat analogous to the energy cascade of turbulence, although with some important differences, already discussed by Hoyle (1953, see also Field et al. 2008), most important of which is the fact that the mass cascade implies that all scales accrete from their parent structures.

In view of the above, the time is ripe for a re-analysis of the evidence and arguments leading to the extended notion that MCs require support, and for a comprehensive formulation of the GHC scenario, which has been presented in scattered form in various papers (Vázquez-Semadeni et al., 2007, 2009, 2010; Vázquez-Semadeni et al., 2017; Heitsch & Hartmann, 2008; Heitsch et al., 2009; Ballesteros-Paredes et al., 2011a, 2015, 2018; Zamora-Avilés, Vázquez-Semadeni & Colín, 2012; Zamora-Avilés & Vázquez-Semadeni, 2014; Kuznetsova et al., 2015, 2017, 2018a, 2018b; Naranjo-Romero et al., 2015; Camacho et al., 2016), including some new calculations, concerning the timescales involved in the process.

The plan of the paper is as follows: In Sec. 2 we start with a critical discussion of the arguments that originally led to the notion that MCs require support (Sec. 2.1), and then continue with a discussion of some problems with the notion of turbulent support (Sec. 2.2). Next, in Sec. 3 we revisit the scenario of hierarchical gravitational fragmentation in an isothermal medium (Hoyle, 1953), and discuss its extension to the case of nonlinear density fluctuations produced by turbulence, in particular the timescales involved in the sequential onset of collapse of different scales. In Sec. 4 we then revisit some particular features and misconceptions about the mechanism of non-homologous gravitational collapse that explain various structural features of observed cores. Next, in Sec. 5, we give a comprehensive description of the GHC model, covering the development of filamentary structure, the increase of the SFR, the formation of massive star-forming regions, the termination of the local SF episodes and the self-regulation of star formation (SF) by feedback, and finally the resulting stellar cluster structure. In Sec. 6 we discuss implications and caveats of the model, and in Sec. 7 we conclude with a summary and some final remarks.

2 Do molecular clouds need to be supported?

In this section we revisit the arguments that originally led to the established notion that MCs need to be supported against gravity, showing that they do not hold in the light of current observational data of the structure and dynamics of MCs, and then we summarize some difficulties with the hypothetical turbulence that would be able to provide support for the clouds.

2.1 Deconstructing the standard objections against global cloud collapse

As mentioned in Sec. 1, one of the first interpretations of the supersonic linewidths observed in MCs (Wilson et al., 1970) was that the clouds are in global gravitational collapse (Goldreich & Kwan, 1974; Liszt et al., 1974). However, this proposal was shortly afterwards rejected by Zuckerman & Palmer (1974, hereafter ZP74) and Zuckerman & Evans (1974, hereafter ZE74), on the basis of two main arguments:

The SFR conundrum: If the clouds were in global collapse, then the star formation rate (SFR) would be much larger than the actual observed value. 2. 2.

The line-shift-absence conundrum: Absorption lines produced by the gas surrounding star-forming regions should be systematically redshifted from emission lines produced at the region itself, since the surrounding gas should be infalling into the star-forming region. Such systematic redshifts, however, were not observed.

In view of these arguments, ZE74 argued that the supersonic linewidths observed in MCs should correspond to small-scale supersonic turbulent motions. Note that the requirment that the motions be small-scale is essential, since it is necessary to avoid cloud-scale systematic motions that would produce the non-detected line shifts between envelope and core regions.

However, the notion of small-scale turbulence (the only one that can resolve the ZE74 conondrums) faces several problems of its own. First, the SFR conundrum has had an alternative solution that dates back to the same times in which it originated, namely that the clouds may be destroyed early in their life cycle by their first stellar products, so that only a small fraction of their mass manages to be converted to stars during the collapse (e.g., Field, 1970; Whitworth, 1979; Elmegreen, 1983; Cox, 1983; Franco, Shore & Tenorio-Tagle, 1994). That is, rather than the clouds being in equilibrium and maintaining a low SFR over their entire lifetimes, the alternative is for SF to constitute a highly intermittent process in space and time, with short and intense bursts of SF that destroy their parent clouds before most of their mass can be converted to stars.

Indeed, recent numerical simulations have shown that the ionizing radiation from massive stars is generally capable of eroding away MCs of masses up to $\sim 10^{5}M_{\odot}$ (e.g., Dale et al., 2012; Colín et al., 2013; Wareing et al., 2017a, b; Haid et al., 2019), although the destruction and dispersal of MCs more massive than that may require additional driving mechanisms, and remains a topic of debate. For example, using a one-dimensional, spherically symmetric model, Murray et al. (2010) proposed that radiation pressure is an essential ingredient for dispersing GMCs in all environments, while protostellar jets may be important only during early evolutionary stages of the clouds, and ionising radiation is effective in Milky-Way-type galaxies, but not in starbursts. Numerical simulations of initially spherical clouds with decaying turbulent velocity fields confirm that photoionising radiation appears incapable of destroying clouds more massive than $10^{5}M_{\odot}$ (e.g., Dale et al., 2012, 2013), although simulations also suggested that radiation pressure would be less effective than simple 1D models predicted, because of the development of Rayleigh-Taylor instabilities that limit the coupling between radiation and matter (e.g., Krumholz & Thompson, 2012, 2013; Reissl et al., 2018). On the other hand, the effectiveness of supernova (SN) explosions in destroying clouds depends strongly on the relative locations of the SNe and the dense gas, although when the SNe are self-consistently positioned, cloud dispersal is generally accomplished for low- to intermediate-mass clouds (e.g., Iffrig & Hennebelle, 2015, 2017; Körtgen et al., 2016; Seifried et al., 2018). Apparently, the destruction of massive clouds requires the combined action of winds and the thermal and radiation pressures (Rahner et al., 2017).

In any case, it is well established observationally that stellar associations and clusters rid themselves of their gas within 5-10 Myr (e.g., Leisawitz et al., 1989; Lada & Lada, 2003; Shimoikura et al., 2018).

On the other hand, the line-shift-absence conundrum of ZP74 is easily explained through geometrical considerations. Essentially, the arguments leading to this conundrum assume that the collapse is roughly spherically symmetric and monolithic, so that the infall motions are coherent, and directed towards a single collapse center, at the geometrical center of the cloud. This assumption is inconsistent with our current understanding of the structure of MCs, which are known to be far from spherically symmetric, and instead consist of an intrincate and inhomogeneous network of filaments and clumps within them (e.g., Bally et al., 1987; Feitzinger et al., 1987; Gutermuth et al., 2008; Myers, 2009; Juvela et al., 2009; André et al., 2010; Henning et al., 2010; Men’shchikov et al., 2010; Molinari et al., 2010; Arzoumanian et al., 2011; Busquet et al., 2013). The central clumps (“hubs”) appear to accrete from the filaments, while in turn the filaments seem to accrete radially from their surroundings (Schneider et al., 2010; Kirk et al., 2013; Peretto et al., 2014; Lu et al., 2018; Williams et al., 2018; Shimajiri et al., 2019; Gong et al., 2018). Thus, the geometry is far from being spherically symmetric, and therefore the accreting gas is not isotropically distributed around the collapse centers (the hubs). In addition, the velocity field is highly complex and chaotic (e.g., Gómez & Vázquez-Semadeni, 2014; Zamora-Avilés et al., 2017; Gómez et al., 2018), so there is no reason to expect a systematic redshift of the absorption lines produced in the gas surrounding the hubs. Instead, the accretion flow is most directly observed as velocity-centroid gradients along the filaments, directed towards the hubs. Indeed, synthetic CO observations of simulations of the regime often show only marginal or no evidence for infall profiles, due to the chaotic motions and perhaps velocity crowding effects, although the line profiles do look similar to observed ones (e.g., Heitsch et al., 2009; Heiner et al., 2015; Clarke et al., 2018). Nevertheless, recent dedicated searches for evidence of infall signatures in CO lines from GMCs have met with success. For example, Schneider et al. (2015) have found the classical combination of self-absorbed and blue-skewed optically thick lines (12CO ($3\rightarrow 2$)) together with centrally-peaked optically thin (13CO ($1\rightarrow 0$)) lines, indicating collapse in the molecular gas surrounding IRDC G28.37+0.07, while Barnes et al. (2018) have measured shifts between the lines of 12CO (tracing gas in the outer parts of the cloud) and 13CO (tracing gas deeper into the cloud) in the pc-scale, massive clumps of the CHaMP survey, finding systematic velocity differentials between the two lines that imply an average mass accretion timescale of $\sim 16$ Myr, consistent with the timescales we discuss in this paper (cf. Sec. 7.1 and Fig. 13).

2.2 Difficulties with turbulent support

It is also very important to note that the ZE74 line-shift-absence conundrum would not only rule out global collapse of clouds, but also MC turbulence in general, as we presently understand it. Indeed, real-world turbulence is characterized by an energy spectrum containing most of the energy at large scales. This implies that the average velocity difference between points in the flow is larger for larger separations between the points, with this trend continuing up to the scale at which the energy is injected (the “energy injection”, or “energy containing” scale; see, e.g., Frisch, 1995).

It is important to emphasize that the larger velocities at large scale refer to coherently-moving structures. That is, they do not only imply a larger velocity dispersion for larger regions; rather, they mainly imply large-scale coherent motions, up to the scale of the clouds themselves. Therefore, the large-scale turbulent motions must shred, compress, or gyrate the clouds, rather than keeping them in a quasi-static state.

Indeed, studies aimed at determining the energy-containing scale in MCs conclude that it corresponds to the largest scales in the clouds (e.g., Brunt et al., 2009), and so the dominant motions in the clouds are those at their largest scales. This implies that the objection to global gravitational contraction by ZE74 would also apply to the large-scale motions that dominate the turbulence in MCs. One concludes that either both global gravitational collapse and large-scale-driven compressible turbulence are ruled out by the ZE74 argument, or else neither one is.

In addition, plain turbulence as we understand it today has problems explaining some observed properties of MCs:

Clouds both with and without SF have similar nonthermal velocity dispersions (e.g., Williams & Blitz, 1998) 2. 2.

Clouds appear to be close to equipartition between their kinetic and gravitational energies, except at low column densities, where large kinetic energy excesses are often found (e.g., Larson, 1981; Heyer et al., 2009; Kauffmann et al., 2013; Leroy et al., 2015). This fact has traditionally been interpreted as virial equilibrium, although Ballesteros-Paredes et al. (2011a) have argued that it rather corresponds to the signature of gravitational collapse at all scales.111In addition, as discussed by Ballesteros-Paredes (2006), virial equilibrium, as it is generally discussed in the literature, assumes that (a) the time derivatives are zero (while in reality, they cannot be measured); (b) clouds have fixed masses and there is neither mass nor energy flux through their boundaries (while in reality, there should be substantial exchange of these quantities between the clouds and their surroundings, first as accretion during cloud growth, and then as return of material during the cloud disruption; see also Ballesteros-Paredes et al., 1999a; Klessen & Hennebelle, 2010; Vázquez-Semadeni et al., 2010; Ibáñez-Mejía et al., 2017); (c) the gravitational potential arises only from the cloud itself and is approximately that of a uniform sphere (while in reality the clouds are far from homogeneous and the environment may dominate the gravitational potential; see also Ballesteros-Paredes et al., 2018); (d) turbulence is continuously driven and strong enough that it can prevent collapse on large scales and promote collapse at smaller scales (while, as we argue in this paper, simulations with self-consistent turbulence driving by the formation of the clouds and by stellar feedback do not support this assumption. Instead, the turbulence generated by the cloud assembly itself cannot support the clouds, while the feedback destroys the clouds rather than supporting them.)

Indeed, the first property rules out driving from internal stellar sources, because this driving cannot explain the turbulence in starless clouds. Conversely, if the driving is external, then strong fine-tuning would be required in order to satisfy the second property, because there is no a priori reason why the turbulence induced in the cloud should acquire precisely an equipartition level. The induced motion might just as well be either too strong, with the effect of dispersing the cloud, or too weak, being insufficient to support the cloud and thus allowing subsequent collapse.

Another view of external driving is that presented by Padoan et al. (2016), who have suggested that MCs in general do not exhibit relations corresponding to near equipartition (of the form of eq. [1]). Instead, they suggest that the clouds exhibit a plain turbulence-scaling relation ${\cal L}\equiv\sigma_{v}/R^{1/2}=$ cst, where ${\cal L}$ is hereinafter referred to as the Larson ratio.222This quantity is sometimes referred to as the velocity scaling or the size-linewidth coefficient.

Padoan et al. (2016) report the scaling ${\cal L}\approx$ cst. for their simulations of the SN-driven ISM, and noted that this is the scaling observed for clouds in the Outer Galaxy Survey (Heyer et al., 2001). They thus argued that this should be considered as the general scaling, because the Outer Galaxy Survey is the largest MC survey available. However, concerns exist about this interpretation. First, the simulation by Padoan et al. (2016) is confined to a (250 pc)3 box, and so the thermal and kinetic energy injected by the SNe cannot escape to high altitudes above the Galactic plane and drive a Galactic fountain, as it should, and so the simulation is most probably overdriven. Second, the Outer Galaxy clouds are likely to have lower column densities and thus to be less influenced by self-gravity than their inner-Galaxy counterparts. This is consistent with the fact that data sets from several nearby galaxies (Leroy et al., 2015) exhibit in general a Heyer-like scaling ${\cal L}\propto\Sigma^{1/2}$, except for the lowest column clouds, which tend to be dominated by turbulence, and have ${\cal L}\sim$ cst., although with a large scatter. This trend is reproduced in simulations of less intense turbulence in the ISM (Camacho et al., 2016; Ibáñez-Mejía et al., 2016).

Note that the scenario of turbulence driven by internal feedback from stellar sources has been assumed by many analytical or semi-analytical models, which have suggested that the SFR can self-regulate and maintain near-virial balance in the clouds (e.g., Norman & Silk, 1980; Franco & Cox, 1983; McKee, 1989; Krumholz et al., 2006; Goldbaum et al., 2011), and by numerical simulations of the effect of feedback from stellar outflows on clump-sized ($\sim 1$ pc) structures (e.g., Nakamura & Li, 2005, 2007; Li & Nakamura, 2006; Carroll et al., 2009, see also the reviews by Vázquez-Semadeni 2011 and Krumholz et al. 2014). However, the full numerical simulations of MCs mentioned above (Dale et al., 2012, 2013; Krumholz & Thompson, 2012, 2013; Colín et al., 2013; Iffrig & Hennebelle, 2015, 2017; Körtgen et al., 2016; Wareing et al., 2017a, b; Seifried et al., 2018), which have included various forms of feedback, generally suggest that the result is to destroy the cloud by a combination of “evaporation” to an ionized phase and dispersal of the remaining cold, dense fragments, rather than maintaining a near-equilibrium state in the cloud. An analytical model that has considered ionizing feedback in non-spherical clouds, with the resulting cloud dispersal, has been presented by Zamora-Avilés, Vázquez-Semadeni & Colín (2012, see also Völschow et al. 2017).

On the other hand, the possibility of external driving has been favored by numerical simulations of cloud formation by converging flows in the warm neutral medium (WNM) (e.g., Vázquez-Semadeni et al., 2007, 2010; Heitsch & Hartmann, 2008; Klessen & Hennebelle, 2010; Micic et al., 2013). However, the turbulence produced in the clouds during their assembly is generally found to be too weak to support the cloud, and global collapse ensues nevertheless. This is essentially a consequence of the facts that a) the turbulence injected by the assembly process remains at a constant level, while the mass of the cloud continues to increase (see also Folini & Walder, 2006), and b) this constant level of the velocity dispersion is only moderately supersonic (sonic Mach numbers ${\cal M}_{\rm s}$ of a few; e.g., Koyama & Inutsuka, 2002; Audit & Hennebelle, 2005; Vázquez-Semadeni et al., 2006, 2007; Banerjee et al., 2009), significantly weaker than the highly supersonic regime (${\cal M}_{\rm s}\ga 10$) known to exist in GMCs (e.g., Heyer & Brunt, 2004).333Note that Klessen & Hennebelle (2010) report that, in their simulations, the internal velocity dispersion of clumps produced in their colliding flow simulations agrees with a Larson scaling, but it should be noted that the vast majority of the clumps they examined had sizes smaller than 1 pc, so they did not probe the GMC regime.

Global hierarchical collapse, on the other hand, naturally explains the properties of the observed motions in the clouds, because the motions are not driven by stellar feedback, and so no stellar population needs to be present in order to attain the observed motions, explaining property (1) above. In addition, the motions necessarily correspond to equipartition with the gravitational energy, since they amount essentially to free-fall, and the free-fall speed is observationally indistinguishable from the virial speed, given typical uncertainties in the measurements (Ballesteros-Paredes et al., 2011a), thus explaining property (2).

It is worth remarking that the ISM and MCs are undeniably highly turbulent, given their very large Reynolds numbers (e.g., Elmegreen & Scalo, 2004). However, the conclusion from this section is that the nonthermal motions observed in MCs are likely to consist of a combination of infall and truly turbulent motions, resulting in what Ballesteros-Paredes et al. (2011a) refer to as “chaotic infall”. The latter may tend to oppose the collapse at large scales (while producing small-scale compressions), but they are likely not as strongly supersonic as they are normally thought to be, since they may actually be seeded by the converging flows that form the clouds, but only to a moderately (rather than strongly) supersonic regime, and then fed by the collapse itself (e.g., Vázquez-Semadeni et al., 1998; Robertson & Goldreich, 2012; Murray & Chang, 2015), therefore being a byproduct of the collapse rather than an external agent capable of preventing it.

3 Hoyle fragmentation revisited

3.1 General considerations and assumptions

If we accept the hypothesis that the supersonic motions observed in the clouds are not dominated by random turbulence (except perhaps those observed in the lowest-column-density clouds; Leroy et al., 2015; Camacho et al., 2016; Traficante et al., 2018a), a simple alternative is to return to the proposal of Goldreich & Kwan (1974), that they are dominated by infall. However, in view of the discussion in Sec. 2.1 above, the infall should not be assumed to be monolithic nor spherically symmetrical. Instead, it can be safely considered to be hierarchical, as we will show in this section. The fundamental underlying mechanism is essentially that outlined by Hoyle (1953), and later laid on a more quantitative mathematical foundation by Hunter (1962, 1964). The mechanism is based on the fact that, for a nearly isothermal medium (and more generally, for any polytropic flow with an effective polytropic exponent $\gamma<4/3$; e.g., Chandrasekhar, 1961; Vázquez-Semadeni et al., 1996), the Jeans mass decreases with increasing density so that, as an isothermal cloud of fixed mass $M$ contracts gravitationally, the number of Jeans masses it contains increases with time, and so, new stages of collapse are initiated as time proceeds, as sequentially smaller masses become unstable. This trend continues until the flow ceases to behave isothermally, as it eventually becomes optically thick, as pointed out by Hoyle (1953) himself, who estimated that the final fragment masses should be $\sim 0.1\,M_{\odot}$, not far from the presently accepted lower limit of the stellar initial mass function.

As appealing as it was, this mechanism later fell in disfavor when it was realized that, for spherical, nearly uniform clouds starting with masses close to the Jeans mass and with linear (small-amplitude) density fluctuations, the largest growth rates correspond to the largest size scales, and so the small-scale perturbations do not have time to complete their collapse before the whole cloud does (Tohline, 1980). However, as mentioned above, the recent numerical simulations of GMC formation including cooling and self-gravity show that fragmentation clearly does occur (e.g., Koyama & Inutsuka, 2002; Audit & Hennebelle, 2005, 2010; Heitsch et al., 2005, 2008a, 2008b; Vázquez-Semadeni et al., 2007, 2010, 2011; Heitsch & Hartmann, 2008; Hennebelle et al., 2008; Banerjee et al., 2009; Colín et al., 2013; Micic et al., 2013; Feng & Krumholz, 2014; Carroll-Nellenback et al., 2014; Fogerty et al., 2016; Körtgen et al., 2016; Wareing et al., 2019).

Two likely reasons why this Tohline barrier, as we will refer to it, is overcome in the numerical simulations are that a) the clouds formed by large-scale coverging flows that coherently trigger phase transitions from the warm to the cold atomic phase quickly acquire many Jeans masses, and b) the assembly mechanism drives moderately supersonic turbulence (${\cal M}_{\rm s}\sim 3$) through the nonlinear thin shell instability (NTSI; Vishniac, 1994), which, aided by the thermal, Kelvin-Helmholz and Rayleigh-Taylor instabilities, causes significantly nonlinear (large-amplitude) density fluctuations (Koyama & Inutsuka, 2002; Heitsch et al., 2006; Vázquez-Semadeni et al., 2006; Banerjee et al., 2009; Wareing et al., 2019). Moreover, the geometry of the flow collisions leads to the formation of flattened rather than spherical clouds. Thus, the main conditions on which the Tohline barrier is based, namely a spherical geometry, linear density fluctuations, and the closeness to the Jeans mass, are violated in the mechanism of GMC formation by converging flows in the warm neutral medium, and thus Hoyle-type gravitational fragmentation appears feasible again.

In the remainder of this section, we first outline the sequential, global hierarchical fragmentation mechanism, and then compute the timescale for the onset of the collapse of objects of different masses. The mechanism proposed can be summarized as follows:

The global, hierarchical gravitational collapse amounts to a mass and energy cascade from large to small scales driven by self-gravity (Field, Blackman & Keto, 2008), so that at each moment in time each structure is accreting from its parent structures and onto its own child substructures. 2. 2.

The nonthermal motions in the gas consist of two main components: a moderately supersonic turbulent background, with typical sonic Mach number ${\cal M}_{\rm s}\sim 3$, and a dominant infall flow. 3. 3.

The global infall flow is hierarchical, and consists of large-scale flows (which are often filamentary; see Sec. 5.3) directed toward the minimum of the large-scale potential wells, on top of which ride flows directed toward local potential minima, therefore causing Hoyle-like fragmentation. This type of flow has been observed numerically (Smith et al., 2011; Gómez & Vázquez-Semadeni, 2014; Zamora-Avilés et al., 2017; Gómez et al., 2018) and suggested observationally (Longmore et al., 2014). The latter authors refer to this type of flow as a “conveyor belt” flow, and we will maintain this nomenclature. 4. 4.

The hierarchical collapse process consists of the sequential destabilization of progressively smaller-mass density fluctuations, as the global collapse causes an increase of the mean density within the cloud, and thus causes the mean Jeans mass in the cloud to decrease. Therefore, progressively smaller-mass turbulent density fluctuations become Jeans-unstable as time proceeds, and begin their own collapse at sequentially later times. 5. 5.

As the cloud’s mean density increases and the average Jeans mass decreases, the cloud’s mass becomes progressively larger than the mean Jeans mass, and the large-scale collapse approaches a pressureless regime. This causes the collapse motions to amplify anisotropies (Lin et al., 1965) and to produce sheets and filaments (Gómez & Vázquez-Semadeni, 2014). 6. 6.

The collapse of each mass scale starts at a finite time during the evolution of the cloud, with a finite initial radius, and from zero local infall speed (Ballesteros-Paredes et al., 2018). Figure 1 schematically illustrates the first step of this hierarchical collapse process. 7. 7.

While the small scales undergo their own collapse, they simultaneously participate of the large-scale collapse; i.e., they fall into the large-scale minimum of the gravitational potential, and so the large-scale collapse amounts to the merger of the locally-collapsing regions.

3.2 Definitions and order-of-magnitude quantities

In this section we will consider a Hoyle-like mechanism assuming spherical geometry which, although higly idealized, illustrates the essential aspects of GHC, and also serves as the basis for further refinements, some of which we will consider in the subsequent sections.

For the sake of calculation, we consider a spherical isothermal cloud that starts out with a fixed mass $M$ just above its Jeans mass. Although, as mentioned in Sec. 3.1, in reality MCs are expected to quickly acquire many Jeans masses, the calculation is made simpler by assuming our cloud contains just over one Jeans mass but can fragment nevertheless and collapses pressurelessly. This assumption also applies to later stages of fragmentation. Generalization to a multi-Jeans mass scenario is trivial. We also neglect the role of accretion onto each level of the hierarchy.

The Jeans mass is in general defined as

[TABLE]

That is, $M_{\rm J}$ is defined as the mass contained in a uniform sphere with density $\rho$, isothermal sound speed $c_{\rm s}$, and radius $L_{\rm J}/2$, where

[TABLE]

is the Jeans length.

The cloud thus starts contracting gravitationally, but, as it contracts and its mean density increases, its mean Jeans mass decreases, so the number of Jeans masses it contains increases with time. Also, because fragments at all levels of the hierarchy are collapsing on their own, after some time (see below), they will all exhibit internal contraction velocities of the order of their free-fall velocity,

[TABLE]

where $M$ and $R$ are the mass and radius of a fragment, respectively, and we have neglected a geometrical coefficient of order unity which, for a uniform-density sphere, is 3/5. However, this velocity is an upper limit, because it assumes that the fragment had an infinite radius when it began contracting. In reality, as stated in item 6 of Sec. 3.1, the local contraction began with a certain finite radius $R_{\rm 0}$ equal to half the local Jeans length. Thus, the actual infall speed can be calculated from energy conservation (Ballesteros-Paredes et al., 2018), writing $E_{\rm k}+E_{\rm g}=E_{\rm tot}$, where $E_{\rm k}=1/2Mv_{\rm g}^{2}$ is the infall kinetic energy, $E_{\rm g}=-GM^{2}/R$ is the instantaneous gravitational energy, and $E_{\rm tot}=-GM^{2}/R_{\rm 0}$ is the gravitational energy of the fragment (neglecting geometrical factors of order unity) at the time when it started its own contraction. Thus, the gravitationally-driven infall speed is given by

[TABLE]

Using the assumptions of sphericity and of a fixed mass $M$, this expression can be written in terms of the column density $\Sigma\equiv M/\pi R^{2}$ and the Larson ratio ${\cal L}$ (cf. Sec. 2.2; Keto & Myers, 1986; Heyer et al., 2009) as

[TABLE]

where the subindex “g” denotes the Larson ratio corresponding to the gravitational velocity given by eq. (5), and $\Sigma_{0}$ denotes the initial column density of the object.

In addition, Ballesteros-Paredes et al. (2018) also considered that the cloud may have started with an initial inertial (i.e., not driven by the self-gravity of the cloud) or “turbulent” velocity which, however, does not necessarily provide support. Instead, it may constitute a large-scale compression driven by a passing supernova shock front, or be a generic compressive velocity field driven a larger-scale agent, such as the gravitational potential of a stellar spiral arm into which the gas may be entering. Ballesteros-Paredes et al. (2018) further assumed that the inertial velocity scales with the size of the clump as

[TABLE]

where, for example, $\eta=1/2$ would correspond to a standard supersonic turbulent scaling, while larger values would allow for the possibility that the inertial motions dissipate in the objects as they contract. As noted by Ballesteros-Paredes et al. (2018), the scaling given by eq. (7) should not be confused with that followed by the component of the velocity driven by self-gravity, which, as pointed out in Ballesteros-Paredes et al. (2011a), is instead expected to follow a relation of the form $\sigma_{v}/R^{1/2}\propto\Sigma^{1/2}$ (cf. Sec. 2.2); i.e., the gravitational velocity depends not only on size, as the turbulent one does, but also on the column density. In this case, adding the gravitational and the inertial components of the velocity in quadrature, the resulting expression for the total Larson ratio is

[TABLE]

Figure 2 shows the evolution of cores following eq. 8, for various cases of the inertial, or “turbulent” velocity in the ${\cal L}$ vs. $\Sigma$ diagram (Keto & Myers, 1986; Heyer et al., 2009), which we refer to as the ${\cal L}$-$\Sigma$ diagram. It is noteworthy that, in all cases, when the collapse is in an advanced stage, the velocity approaches the limit given by the free-fall speed, eq. (4).

It is also convenient to write the virial paramenter, $\alpha_{\rm vir}\equiv 2E_{\rm k,tot}/|E_{\rm g}|$ (Bertoldi, & McKee, 1992), in the case that the kinetic energy $E_{\rm k,tot}$ contains a contribution from truly turbulent random motions, with velocity dispersion $\sigma_{\rm turb}$, and another from gravitationally-driven infall motions, with characteristic velocity $v_{\rm g}$, given by eq. (5). Noting that, in the definition of $\alpha_{\rm vir}$, $E_{\rm g}=1/2Mv_{\rm ff}^{2}$, and using eqs. (4) and (5), we obtain

[TABLE]

We will use this expression in our discussion of the inverse relation between $\alpha_{\rm vir}$ and mass in Sec. 6.4.1.

To conclude this section, it is important to note that, in this basic picture, we are neglecting all geometrical aspects, such as the fact that, in realistic numerical simulations of MCs, the collapse is systematically observed to develop filamentary structures (Gómez & Vázquez-Semadeni, 2014), and the collapse of filaments involves longer timescales than the free-fall time for spherical symmetry (Toalá et al., 2012; Pon et al., 2012).

3.3 The sequential onset of collapse of progressively lower-mass

objects

We now calculate the timescales for sequential destabilization of progressively smaller masses of the typical (rms) density fluctuation. Note that below we also calculate a different timescale: that of the most extreme density fluctuations, which are the first to occur.

The calculation of the destabilization timescale for the typical density fluctuation requires a specific model for the nature of the density fluctuations, which, as stated in Sec. 3.1, need to be either nonlinear, and/or contain many Jeans masses in order to be able to continue fragmenting. Here we perform a proof-of-concept calculation assuming that the fluctuations are generated by moderately supersonic turbulence, which, as stated in assumption 2 above, we generically describe by a typical Mach number ${\cal M}_{\rm s}\sim 3$, so that the typical turbulent density fluctuations have an rms amplitude of order $\rho_{{\rm rms}}/\rho_{0}\sim{\cal M}_{\rm s}^{2}\sim 10$, where $\rho_{0}$ is the mean density in the medium.

Let $M_{0}$ be the total mass of the cloud. The time for an arbitrary fluctuation of mass $M<M_{0}$ to become unstable is the time for the Jeans mass in the typical (rms) density fluctuation to become equal to its actual mass $M$; that is, for $M_{\rm J,rms}(t)=M$. We now proceed to calculate this time. The evolution of the radius as a function of time during free-fall is a transcendental equation. Therefore, in what follows, it will be convenient to use the analytic fit to this evolution provided by Girichidis et al. (2014),

[TABLE]

where $\tau\equiv t/t_{{\rm ff},0}$ is the time in units of the free-fall time for the initial density ($\rho_{0}$) of the object, $t_{{\rm ff},0}$, $R_{0}$ is the initial radius, and $a\approx 1.8614$ is a parameter for which the fit remains within 0.5% of the actual solution during the whole collapse time.

A very important feature of the free-fall process is that it is extremely nonlinear, proceeding very slowly at the beginning, and accelerating enormously towards the end. Girichidis et al. (2014) note that, after 50% of the free-fall time, the radius has only decreased by $\sim 16\%$, and that after 99% of the collapse time, the radius has only decreased by one order of magnitude.

Given expression (10) for the temporal evolution of the radius, we can write an expression for the evolution of the mean density,

[TABLE]

and, from this, an equation for the Jeans mass at the instantaneous mean density; that is, the instantaneous mean Jeans mass in the cloud:

[TABLE]

Now, as stated above, the typical rms density fluctuation is given by

[TABLE]

and the Jeans mass in this fluctuation is

[TABLE]

This equation can be then inverted to find the time required for a mass $M=M_{\rm J,rms}$ to become unstable. We find for this time, in units of the initial free-fall time in the cloud,

[TABLE]

where $\mu\equiv{\cal M}_{\rm s}M/M_{\rm J}(\rho_{0})$ is the fragment mass normalized to the initial fluctuation Jeans mass in the cloud. Figure 3 shows the dependence of this destabilization time (in units of the initial free-fall time) for mass $M$ (in units of the initial Jeans mass) for various values of the Mach number. We consider only masses $\leq M_{\rm J,0}/2$ assuming that fragmentation starts when there are at least two Jeans masses in the cloud.

We can see that, for larger Mach numbers, smaller masses can begin to collapse during the collapse of the largest scale.

It is important to note that the above calculation is only an approximation intended as a proof-of-concept for the sequential onset of collapse of the hierarchy of mass scales in a cloud, and probably provides only an upper limit to this time. A more precise calculation should consider several additional factors that can shorten the time for the onset of collapse of structures of mass $M$, such as:

• •

This calculation refers only to the first stage of fragmentation for any given mass. We have neglected the subsequent fragmentation of a fragment.

• •

The collapse of a structure embedded in an also-contracting background takes a shorter time than the standard collapse time for isolated objects (Toalá et al., 2015).

• •

In reality, there is a distribution of density fluctuations, and therefore fluctuations of a higher density are expected to exist (if albeit rare), so that their free-fall times are shorter.

Nevertheless, the main point of the above calculation is to show the fact that progressively smaller masses become unstable at progressively later times.

3.4 The time to reach the $n$-th level of fragmentation

Another relevant quantity that can be computed for the fragmentation process is the time to reach the $n$-th level of fragmentation. This can be done by assuming that a new fragmentation stage begins every time an object undergoing gravitational contraction in the previous stage attains a number $n_{\rm MJ}$ of local Jeans masses as its density increases and the local Jeans mass decreases. This is consistent with the empirical finding by Guszejnov et al. (2018) in their numerical simulations that a structure must contain at least 3 Jeans masses in order for it to fragment instead of collapsing monolithically. Thus, the Jeans masses at stages $n$ and $n-1$ of the hierarchy are related by

[TABLE]

As in the previous section, we consider that the fragmentation hierarchy starts at the Jeans mass in the rms turbulent density fluctuation which, according to eq. (14), at the onset of global collapse, is given by $M_{\rm J,rms}(\tau=0)=M_{\rm J,0}/{\cal M}_{\rm s}$. We can then write the successive steps of the hierarchy as

Level 1:

[TABLE]

Level 2:

[TABLE]

Level $n$:

[TABLE]

On the other hand, if $\tau_{n}$ is the time at which the $n$-th fragmentation level is reached, then, from eq. (14) we can write

[TABLE]

Therefore, solving for $\tau_{n}$ from equations (17) and (18), we find

[TABLE]

which gives the necessary time to reach the $n$-th level of fragmentation. The left panel of Fig. 4 shows plots of this expression for $n_{\rm MJ}=2$, 3, 4, and 5 (solid, dashed, dot-dashed, and dash-dot-dotted lines, respectively). The right panel of this figure shows the reverse plot, of the number of fragmentation levels reached at a given elapsed time.

It is noteworthy that this result does not depend on the Mach number of the turbulence. We discuss insights and implications of this result in Sec. 6.6.1.

3.5 The time for the onset of the first fragment contraction

Until now, we have considered the times for the onset of gravitational contraction for the typical density fluctuations, of amplitude proportional to ${\cal M}_{\rm s}^{2}$. However, the first fragments that will become Jeans unstable are not those of typical amplitude, but rather the most extreme ones, as they have the shortest free-fall times. We can compute the time for the onset of the first collapses by noticing that they can only occur when the total mass at the highest densities $M(n\geq n_{\rm SF})$ is at least equal or higher than the Jeans mass at those densities $M_{\rm J}(n_{\rm SF})$. That is, we write the condition for the occurrence of the first collapse as

[TABLE]

where we have denoted by $n_{\rm SF}$ the density for which the free-fall time is so much shorter than the average free-fall time in the cloud that it can be considered instantaneous. That is, once condition (20) is met, the collapse of this region will occur essentially immediately.

We assume a lognormal density PDF, appropriate for supersonic, isothermal turbulence (Vázquez-Semadeni, 1994; Padoan et al., 1997; Passot & Vázquez-Semadeni, 1998; Federrath, Klessen & Schmidt, 2008), of the form

[TABLE]

where $s\equiv\ln(\rho/\bar{\rho})$, $b$ is the compressibility parameter (Federrath, Klessen & Schmidt, 2008), $\sigma_{s}^{2}=\ln(1+b^{2}{\cal M}_{\rm s}^{2})$, and $\langle s\rangle=-\sigma_{s}^{2}/2$. In this case, the mass $M(n\geq n_{\rm SF})$ is given by (Krumholz & McKee, 2005)

[TABLE]

Note that, in fact, this mass estimate is strictly a lower limit, since the PDF contains no spatial information, and it is thus possible that not all the mass above $n_{\rm SF}$ lies in a single connected object in the medium. However, since we are considering the global density maximum of the flow, it is reasonable to assume that the next densest fluid parcels belong to the same object. With this caveat in mind, the left panel of Fig. 5 shows plots of the ratio $M(n\geq n_{\rm SF})/M_{\rm J}(n_{\rm SF})$ at a fixed value of $M_{\rm tot}=10^{4}M_{\odot}$, for various values of the rms Mach number, assuming $\langle n\rangle=10^{3}\,{\rm~{}cm}^{-3}$ and $n_{\rm SF}=10^{6}\,{\rm~{}cm}^{-3}$, and the right panel shows this ratio for various values of $M_{\rm tot}$ assuming ${\cal M}_{\rm s}=4$, and the same values of $\langle n\rangle$ and $n_{\rm SF}$. The horizontal line shows the value of $M_{\rm J}(n_{\rm SF})$. The time of first collapse for these various cases is given by the intersection of the curves and the horizontal identity line. It can be seen that, in general, the time of first collapse occurs relatively late in the collapse at the mean density of the cloud,, at $t\sim 0.65$–$0.9\,t_{{\rm ff},0}$, where we recall that $t_{{\rm ff},0}$ is the free-fall time at the initial mean density. This is qualitatively consistent with the evolution of numerical simulations, where the first collapsed objects appear several megayears after the onset of the large-scale contraction, but a few megayears before it terminates (e.g., Vázquez-Semadeni et al., 2007, 2010; Heitsch & Hartmann, 2008; Carroll-Nellenback et al., 2014).

4 The nature of the gravitational contraction

In this section we now discuss a few properties of non-homologous gravitational contraction that are often overlooked.

4.1 The non-homologous, non-simultaneous nature of the collapse

A frequent assumption about the process of gravitational contraction is that it proceeds homologously; i.e., that it proceeds maintaining uniform (although not constant in time) density (e.g., Chandrasekhar, 1939). This mode of collapse is accomplished if the infall velocity increases linearly with radius, so that all the shells of a uniform-density spherical structure reach the center at the same time. This idealized solution, however, never occurs in practice in media containing multiple Jeans masses, because any density perturbations present have shorter free-fall times, and thus collapse faster, than the background medium. The requirement of having multiple Jeans masses allows for there to exist a sufficiently large range of unstable scales (larger than the Jeans length) that the reduction of the growth rate at scales just over the Jeans length (cf. Sec. 2.1) is negligible. This can be illustrated by recalling the dispersion relation for the Jeans instability analysis which, for linear perturbations of the form $\rho_{1}(x,t)=\epsilon\rho_{0}\exp\left[-i\left(kx-\omega t\right)\right]$, reads (e.g., Shu, 1992)

[TABLE]

where $\omega$ is the growth rate, $k$ is the wavenumber, and $\rho_{0}$ is the mean density of the medium. As is well known,

[TABLE]

is the wavenumber corresponding to the Jeans length, for which the growth rate is zero; i.e., the restoring force from the pressure gradient balances the gravitational pull in the marginal equilibrium case. At the other extreme, the growth rate is maximum for $k\rightarrow 0$; i.e., at the largest physical scales, at which the growth rate approaches the free-fall value, $\omega\rightarrow\sqrt{4\pi G\rho_{0}}$, so that, in the linear regime, all large-scale perturbations grow at nearly the same rate.

This linear argument, however, neglects the nonlinear effect of the reduction of the free-fall time in the perturbation due to the locally enhanced density. Therefore, when the dynamic range of the large-scale perturbations (for which the pressure gradient is negligible) is sufficiently large, the reduction of the free-fall time at the local density perturbations allows them to outgrow the background. Indeed, it has been recently shown numerically by Guszejnov et al. (2018) that, for numerical boxes containing more than $\sim 3$ Jeans masses, and in the presence of nonlinear turbulent density fluctuations, fragmentation does occur, while, when the numerical box containing less than $\sim 3$ Jeans masses, a single focused, monolithic collapse mode ensues. Moreover, they showed that the number of fragments is independent of the Mach number of the turbulence, which suggests that the fragmentation is of purely gravitational origin.

The above discussion suggests then that the collapse of multi-Jeans-mass clouds in the presence of nonlinear density fluctuations is highly non-homologous, and thus the density contrast in the cloud tends to increase over time, because the denser parts contract faster than the lower-density ones, and quickly overtake them.

An alternative way of seeing this is by noting that the standard similarity treatment of gravitational contraction (e.g., Larson, 1969; Penston, 1969; Shu, 1977; Hunter, 1977; Whitworth & Summers, 1985) proceeds by defining a similarity variable $\xi\equiv r/c_{\rm s}t$, where $r$ is the radial coordinate, $t$ is the time, and $c_{\rm s}$ is the sound speed. The essence of such an analysis is that one radial position in a collapsing sphere at a given time is equivalent to a different radial position at some other time. This can be interpreted as meaning that, at a given time $t$, the gas at a higher density nearer to the center is at a more advanced evolutionary stage (in its approach to the collapse center) than the gas at lower density further out.

In standard spatial and temporal coordinates, this implies that the gas is flowing from low- to high-density regions, so that, as the gas density increases, a gas parcel moves from a more distant to a more nearby location around the collapsing center. Moreover, this process takes a finite amount of time.

The non-homologous nature of the collapse has various implications, most important of which is that, contrary to the implicit assumption in the idealized reasoning of Zuckerman & Palmer (1974), not all gas in a cloud reaches the singularities (i.e., the stars) at the same time. Instead, the densest parts of the cloud terminate their collapse earlier. This is at the foundation of models for the SFR in clouds, that assume the instantaneous value of the SFR in a turbulent cloud, characterized by a certain density probability distribution (PDF), is given by the mass in fluctuations with densities above a certain threshold, divided by a time of the order of the free-fall time for the threshold (e.g., Krumholz & McKee, 2005; Padoan & Nordlund, 2011; Hennebelle & Chabrier, 2011; Federrath & Klessen, 2012; Zamora-Avilés, Vázquez-Semadeni & Colín, 2012; Völschow et al., 2017; Burkhart, 2018), so that the stars formed by the earliest collapses may prevent the collapse of the rest of the material. We return to the application of this reasoning to the GHC scenario in Sec. 5.4.

4.2 The difference between the pre- and protostellar collapse

regimes

Another feature of the gravitational collapse process is that the pre- and protostellar stages have important qualitative differences in the shape of the radial infall velocity profile. As it is well known from isothermal, spherical collapse calculations (e.g., Larson, 1969; Penston, 1969; Hunter, 1977; Whitworth & Summers, 1985; Naranjo-Romero et al., 2015), the collapse process consists of two stages, separated by the formation of a singularity (i.e., a protostar). These correspond to the “prestellar” and the “protostellar” stages of cores.

During the prestellar stage, the gravitational contraction of isothermal, spherically symmetric objects with a central density enhancement asymptotically approaches a radial configuration consisting of two distinct regions (see, e.g., Fig. C1 of Larson (1969) or panels (a)-(c) in Fig. 5 of Whitworth & Summers (1985)):

Prestellar stage:

[TABLE]

[TABLE]

Note that the infall speeds are negative because they point in the negative radial directions. Note also that the transition between the “inner core” and the “envelope” occurs roughly at a radius equal to the Jeans length corresponding to the central density and temperature (Naranjo-Romero et al., 2015; Keto et al., 2015).

The transition between the two regions occurs around a radius equal to the Jeans length at the central density (e.g., Keto & Caselli, 2010), and approaches the core’s center. The uniform central density inside this radius increases without limit with time, becoming infinite at the time of protostar (singularity) formation, at which time the transition radius becomes zero. That is, at the time of singularity formation, the density profile has an $r^{-2}$ shape at all radii, resembling a singular isothermal sphere, but with nonzero, uniform, and supersonic infall speed at all radii (Naranjo-Romero et al., 2015).

Conversely, during the protostellar stage, the density profile increases monotonically inwards at all radii, diverging toward the core’s center, where a finite-mass object (the protostar) is located. Two radial regions still exist in the core, but they are now mediated by a shock front, located at a radius $r_{\rm s}=c_{\rm s}t$, where $t$ is the time since the formation of the singularity. The density and infall velocity profiles in the two radial regions in the protostellar stage have asymptotic regimes given by

Protostellar stage:

[TABLE]

[TABLE]

It is important to notice that this solution differs strongly from the canonical “inside-out collapse” (IOC) (Shu, 1977). The IOC solution results from assuming that the initial condition for the collapse is a hydrostatic singular isothermal sphere (SIS). This assumption requires the pre-stellar contraction to be quasi-static, with the clump’s weight being balanced by some other force(s), such as magnetic forces or turbulent pressure, and a tiny bit of contraction allowed by, for example, ambipolar diffusion (AD) in the case of magnetic support. Although Shu (1977) gave a number of arguments in favor of this prestellar quasistatic contraction, Hunter (1977) pointed out that numerical simulations in general tend to be more closely described by Larson-Penston solution than by Shu’s solution. This can be understood because, in practice, an SIS configuration is unattainable, since it is the most unstable possible hydrostatic equilibrium in spherical symmetry, as follows from the fact that it is the limit of unstable Bonnor-Ebert (BE) spheres when the ratio of the central-to-peripheral density tends to infinity. As discussed by Whitworth et al. (1996) and Vázquez-Semadeni et al. (2005), such an unstable equilibrium is not expected to arise spontaneously within a significantly turbulent molecular cloud. Moreover, magnetic support is also not expected to provide sufficient support, as GMCs are now known to be magnetically supercritical in general (Crutcher, 2012), and so their gravitational contraction must proceed dynamically rather than quasistatically. Finally, turbulent support that dissipates gradually, as often assumed (e.g., Goodman et al., 1998; Bergin & Tafalla, 2007; Pineda et al., 2010), is not feasible because at every scale, the dominant modes of the turbulence are of the same scale, and thus their effect is not to maintain the structure in near-equilibrium, but rather to shear, compress or expand the structure, as discussed in Sec. 2.2.

A common argument against the possibility of cores being in the process of dynamical contraction already during their prestellar stage is that the final parts of this stage produce supersonic infall speeds at the core’s envelope, while typical estimators of the infall speed, such as blue-excess line profiles of moderately optically thick transitions often imply subsonic infall velocities (e.g., Lee et al., 2001; Campbell et al., 2016). However, as shown by Loughnane et al. (2018), the fact that the largest infall speeds in collapsing cores occur at radii where the density is already decreasing (see eq. 24) may cause a systematic underestimation of the infall speed by these indicators.

5 The complete scenario

In this section we now compile all the evolutionary stages and byproducts of the GHC scenario, as they have been reported in several papers over the last decade.

5.1 Cloud formation

Observationally, “clouds” are often defined as regions detected in some tracer, such as 12CO. However, in numerical and analytical studies, clouds are often defined as cold density enhancements over some warm, diffuse background, which may constitute a different thermodynamic phase of the gas. Here we will adopt this definition, so that clouds can contain both cold atomic and molecular gas.

Under solar galactocentric radius conditions, the azimuthally- and vertically-averaged molecular mass fraction is only 10-20%, and most of the Galactic disk volume is occupied by warm atomic gas, except in the spiral arms, where the gas is mostly molecular (Koda et al., 2016). At the turn of the century, several studies showed that dynamic compressions in the warm atomic phase can nonlinearly trigger a phase transition to the cold atomic phase (e.g., Hennebelle & Pérault, 1999; Koyama & Inutsuka, 2002), while generating turbulence in the cold, dense layer that forms as a result of the compression by the combined action of the nonlinear thin-shell (Vishniac, 1994), thermal (Field, 1965) and Kelvin-Helmholtz instabilities (see Heitsch et al., 2006, for a summary of the instabilities active in cloud formation). This turbulence consists mostly of the velocity dispersion of the cold cloudlets formed by the thermal instablity (TI), which is a fraction of the sound speed in the warm phase. Since the latter is $\sim 10\times$ the sound speed in the cold phase, the cloudlets move with a velocity dispersion that is moderately supersonic (Mach numbers of a few) with respect to the sound speed in the cold phase (e.g., Koyama & Inutsuka, 2002; Vázquez-Semadeni et al., 2006; Hennebelle & Audit, 2007; Banerjee et al., 2009). In particular, Banerjee et al. (2009) showed that, if the conditions of thermal- and ram-pressure balance are imposed separately between the warm and cold phases, then the turbulent Mach numbers in the warm and the cold phases are similar.

Finally, an important result is that, at the earliest stages, the condensed layer is very thin, with column densities $\sim 10^{19}{\rm~{}cm}^{-2}$, and therefore must consist almost exclusively of cold atomic gas, thus being highly consistent with the conclusion of Heiles & Troland (2003) that CNM clouds appear to be large, thin sheets, with aspect ratios of up to $\sim 100$ (Vázquez-Semadeni et al., 2006). Thus, GMCs may start their lives as thin sheets of cold atomic gas.

It is important to remark that this formation mechanism implies that, even though the clouds do constitute a cold, dense phase in the ISM and are surrounded by a warm, diffuse substrate, there is a continuous process of accretion of gas from the warm phase into the clouds. The gas suffers a phase transition from the warm to the cold phase upon entering the clouds, and thus the cloud boundaries are phase transition fronts (Banerjee et al., 2009), but no restriction to the flow of gas across these fronts exists (Ballesteros-Paredes et al., 1999a; Hennebelle & Pérault, 1999; Koyama & Inutsuka, 2002; Audit & Hennebelle, 2005; Heitsch et al., 2005; Vázquez-Semadeni et al., 2006; Banerjee et al., 2009). This kind of boundary is in sharp contrast with the classical picture that the cloud boundaries are contact discontinuities, where the warm phase “confines” the cold one, but there is no motion across the boundaries. This also implies that, in general, cloud masses vary over time rather than being fixed. During the early epochs of the clouds, their masses in general grow over time, until they begin to be eroded by their newly formed stars (Sec. 5.4). Observational evidence for this mass evolution of the clouds has been presented by Kawamura et al. (2009), who showed that GMCs in the Large Magellanic Cloud with more evolved stellar populations tend to have larger masses.

5.2 Onset of global gravitational contraction

An often neglected fact is that, as pointed out by Gómez & Vázquez-Semadeni (2014), when the atomic gas undergoes a transition from the warm, diffuse phase to the dense, cold one, the temperature decreases by a factor $\sim 100$ and the density increases by the same factor, so that the Jeans mass, $M_{\rm J}\propto\rho^{-1/2}T^{3/2}$, decreases by a factor $\sim 10^{4}$, from $M_{\rm J}\approx 7\times 10^{7}\,M_{\odot}$ at $T=7000$ K and $n=0.3{\rm~{}cm}^{-3}$, to $M_{\rm J}\approx 7\times 10^{3}\,M_{\odot}$ at $T=70$ K and $n=30{\rm~{}cm}^{-3}$.

Another phase transition occurs when the gas becomes molecular (e.g., Koyama & Inutsuka, 2000), which, at the same thermal pressure, typically has $T=10$ K and $n=200\,{\rm~{}cm}^{-3}$, and furthermore the number density $n$ drops by a factor of 2, so that the Jeans mass is another two orders of magnitude lower, at $M_{\rm J}\approx 30\,M_{\odot}$, in this gas. That is, in molecular gas the Jeans mass is six orders of magnitude lower than in the WNM. This implies that the phase transition from WNM to molecular gas provides a mechanism for creating strongly Jeans unstable gas, if the transition is coordinated by large-scale compressions that coherently produce large masses of cold, dense gas. As mentioned in Sec. 5.1, the initial step of condensation produces thin sheets of CNM which in the simulations are observed to fragment into a network of clumps and filaments (see, e.g., Fig. 10 of Vázquez-Semadeni et al., 2006). It is this ensemble of clumps and filaments that is observed to be globally gravitationally unstable in numerical simulations of cloud formation and collapse (Vázquez-Semadeni et al., 2007, 2010; Heitsch & Hartmann, 2008; Colín et al., 2013; Carroll-Nellenback et al., 2014). In addition, those simulations also show that, in general, the turbulence (mostly clump-to-clump velocity dispersion) generated in the clump ensemble (“the cloud”) produced by the large-scale compressions is not sufficient to support it against its own self-gravity once it becomes Jeans unstable, and therefore the cloud begins to undergo gravitational contraction.

5.3 Molecule and filament formation

It has been known for a few decades (Franco & Cox, 1983; Hartmann, Ballesteros-Paredes, & Bergin, 2001; Bergin et al., 2004) that the column density in the cold phase at which the gas becomes dominated by gravity is very similar (again, for typical solar neighborhood conditions) to that required for H2-molecule self-shielding. Thus, as the gas begins to contract gravitationally, as described in the previous section, it also begins to become molecular, at column densities $N\sim 10^{21}\,{\rm~{}cm}^{-2}$. It is important to recall that the presence of molecular gas does not appear necessary to allow the gas to cool and thus keep contracting; at densities lower than $10^{4}{\rm~{}cm}^{-3}$, C+ fine-structure emission appears to provide the necessary cooling (Glover & Clark, 2012).

Thus, when the gas begins to contract, it is likely to consist mostly of cold atomic gas, and is expected to gradually become molecular during the contraction process. In turn, the transition to molecular contributes to the second reduction stage of the Jeans mass discussed in Sec. 5.2, and thus strengthens the gravitational contraction by further reducing the temperature (see the phase diagrams in Fig. 8 of Koyama & Inutsuka, 2000). The late appearance of molecular gas during the collapse has been observed in numerical simulations of cloud contraction including a prescription for the formation of H2 and 12 CO molecules (e.g., Heitsch & Hartmann, 2008; Heiner et al., 2015). Thus, a “molecular cloud” does not form suddenly, but rather gradually, so that, as the cold gas cloud contracts, its molecular content increases over time, as described by Vázquez-Semadeni et al. (2018). Nevertheless, in the simulations, the gravitational potential driving the collapse is that of the combined molecular and cold atomic gas. This fact must be taken into account when the gravitational binding of clouds is considered.

In addition, as pointed out by Gómez & Vázquez-Semadeni (2014), the fact that the Jeans mass decreases continuously during the gravitational contraction as long as the gas remains nearly isothermal, implies that the thermal pressure becomes progressively less important, and the collapse proceeds in an increasingly pressureless manner. In turn, this implies that the collapse tends to enhance any anisotropies initially present in the cloud, since the collapse proceeds faster along the shortest dimensions of the cloud, so that triaxial ellipsoids evolve into sheets, and sheets into filaments (Lin et al., 1965). Moreover, since the cloud is likely to have started out with a sheet-like geometry due to its formation at the interface of converging flows, the pressureless collapse is expected to produce filamentary structures, as verified in numerical simulations of collapsing sheet-like clouds (e.g. Burkert & Hartmann, 2004; Hartmann & Burkert, 2007; Vázquez-Semadeni et al., 2007; Heitsch et al., 2008b, 2009; Gómez & Vázquez-Semadeni, 2014; Carroll-Nellenback et al., 2014).

In particular, Gómez & Vázquez-Semadeni (2014, see also Smith et al. 2011) observed the formation of filaments in numerical simulations of contracting clouds, and pointed out that the filaments constitute the collapse flow itself from the large to the small scales, in which the velocity field is such that material from the sheet-like cloud flows onto the filaments roughly perpendicularly to them. Then the flow changes direction at the filament, becoming longitudinal and directed towards cores in the filaments, or at filament intersections. That is, the filaments are akin to rivers, funneling the gas from the cloud to the dense cores where stars form, as illustrated in the left panel of Fig. 6, where it is seen that the velocity field around the filament is directed mainly toward the filament, while inside the filament the flow is mostly longitudinal, pointing toward the dense core.

It is important to note that the change in direction of the flow along the radial direction occurs smoothly, without the formation of noticeable shocks. This can be observed in the right panel of Fig. 6, which shows that the inwards radial velocity decreases smoothly towards the filament’s axis, while the longitudinal component grows also smoothly. No jump in the radial velocity is observed at the axis. This is similar to the smooth nature of the velocity field in a collapsing core during its prestellar stage, since a shock does not form until the time of the formation of the singularity (i.e., the protostar; e.g., Whitworth & Summers, 1985). In the filament, no singularity forms precisely because the flow changes direction as it enters the filament and is redirected towards the core. Therefore, the filament is “drained” by the longitudinal flow, and in fact, the simulations indicate that it does not exceed a certain maximum central density, which probably depends on the total amount of mass in the contracting potential well. Thus, in a sense, the filament becomes “frozen” in a state analogous to the prestellar stages of core collapse, albeit in cylindrical geometry.

5.3.1 The magnetic field in the filaments

An implication of the mechanism of filament formation by anisotropic gravitational contraction has been investigated by Gómez et al. (2018), who found, in magnetized simulations, that the accretion flow from the cloud onto the filaments drags the magnetic field lines along the flow, causing it to be oriented mostly perpendicularly to the filaments around them, as observed (e.g., Palmeirim et al., 2013; Cox et al., 2016). However, the longtudinal flow along the filaments onto the main hubs bends the field inside the filament, until the ram pressure of the longitudinal flow is balanced by magnetic diffusion, causing a characteristic “U” shape. The structure seen in the simulations is qualitatively consistent with that observed in polarization maps (e.g., Planck Collaboration et al., 2016).

5.3.2 A quasi-stationary state in the filaments?

In fact, an argument can be made that gravitationally-formed filaments should approach a quasi-stationary state in which the mass flux from the cloud onto the filament is balanced by the flux from the filament to the core, which should vary only on the characteristic variability timescale for the accretion flow onto the filament. This can be seen as follows. Let us assume that, at some time, the filament has a certain charateristic radius $R$, longitudinal velocity $v_{\parallel}$, and mean density $\bar{\rho}_{\rm f}$. The longitudinal mass flux is then $F_{\|}=\pi R^{2}v_{\parallel}\bar{\rho}_{\rm f}$. Let us now assume that the perpendicular flux, $F_{\bot}$, exceeds $F_{\|}$. This will cause the filament to accumulate mass, and thus cause an increase in its mean density and/or radius. Therefore, $F_{\|}$ must increase as well, provided the longitudinal velocity does not decrease. Conversely, let us now assume that $F_{\bot}<F_{\|}$. In this case, the filament loses mass, causing its density and/or radius to decrease, and so $F_{\|}$ decreases. In both cases, then, the change in $F_{\|}$ makes it approach $F_{\bot}$.

This argument suggests that the condition $F_{\|}=F_{\bot}$ may be an attractor for the filament, which may then approach a quasi-stationary flow regime. Note that we refer to it as a quasi-stationary regime, because we have implicitly assumed that $F_{\|}$ has time to adjust to the value of $F_{\bot}$; i.e., that $F_{\bot}$ varies slowly in comparison to the variation timescale of $F_{\|}$. This is a reasonable assumption, since $F_{\bot}$ is the result of the large-scale accretion flow, while $F_{\|}$ represents a smaller-scale accretion flow. But nevertheless, the quasi-stationary regime must vary on the variability timescale of $F_{\bot}$. Numerical experiments exploring the approach to this quasi-stationary regime will be presented elsewhere (Naranjo-Romero et al., in prep.).

If the central density of the filament saturates because of the establishment of this quasi-stationary regime, then no strong shocks arise in the central parts of the filament due to the accretion flow (except at the cores, which represent local collapse sites in supercritical filaments). The absence of a strong central shock should thus be a unique signature of filaments formed by gravitational contraction, in contrast to filaments formed by strongly supersonic turbulence. This may be used as a diagnostic to distinguish the origin of MC filaments. The only shocks present in the filaments may be weak shocks due to the residual turbulence in the filamentary flow, which is likely to be only mildly supersonic at most in MC filaments.

5.4 Onset and acceleration of the star formation activity

5.4.1 Can turbulence alone induce local collapses?

Numerical simulations of the contraction of MCs after their formation as cold atomic clouds including sink particles (Vázquez-Semadeni et al., 2007, 2010; Colín et al., 2013; Carroll-Nellenback et al., 2014; Feng & Krumholz, 2014) systematically show that the sink particles appear late in the evolution of the cloud, typically $\sim 5$–10 Myr (a few free-fall times, due to geometrical effects; Toalá et al., 2012; Pon et al., 2012) after the onset of global contraction in the cloud, which in turn occurs $\sim 10$ Myr after the first appearance of cold atomic gas triggered by compressions in the warm atomic gas.444 Note that newly formed cold atomic clouds in simulations with physical conditions characteristic of the Solar neighborhood, typically take up to 10 Myr to become Jeans unstable. Note also that the shorter timescales of 1–3 Myr for collapse calculated by Hartmann, Ballesteros-Paredes, & Bergin (2001) correspond to the duration of the collapse process for clouds that are already molecular and Jeans unstable. Instead, the timescales referred to here correspond first to the timescale for cold atomic gas to become Jeans unstable, and then to the duration of the collapse starting out from cold atomic gas conditions, which is typically 10$\times$ warmer and less dense than the molecular gas considered by Hartmann, Ballesteros-Paredes, & Bergin (2001). However, the appearance of these first stars occurs $\sim 5$ Myr before the large-scale collapse has advanced sufficiently to form massive stars.

This clearly implies that the primordial turbulence injected to the clouds by various instabilities during their assembly (Sec. 5) is too weak to directly generate density fluctuations of large enough amplitude for the local Jeans mass to become smaller than the fluctuation mass, triggering local collapse. This is true even in simulations starting with turbulent fluctuations with kinetic energies comparable or even larger than the gravitational energy (e.g., Bate et al., 2003; Bate, 2009b; Dale et al., 2012; Howard et al., 2018; Grudić et al., 2018). These simulations also do not form stars until the turbulence has been dissipated, and gravitational contraction has set in. That is, the initial turbulence does not form stars by itself. In fact, Clark & Bonnell (2005) have suggested that isothermal turbulent compressions do not effectively reduce the local Jeans mass, because they are preferentially one-dimensional.

Instead, as described in Secs. 3.3 and 3.5, the turbulent density fluctuations begin to become unstable at later times, when the average Jeans mass in the cloud has decreased sufficiently due to the large-scale collapse, so that the masses of the local fluctuations can match it, at which point local collapses start. However, as also described in those sections (see also Vázquez-Semadeni et al., 2009), the turbulent density fluctuations are able to terminate their collapse earlier than the whole cloud because they have nonlinear amplitudes, and therefore they have significantly shorter free-fall times than that of the whole cloud.

5.4.2 Acceleration of the star formation activity

As the average Jeans mass in the parent cloud decreases, the fraction of mass at high enough densities increases, implying an acceleration of the SF process; i.e., an increase in the SFR. This acceleration of the SF continues until a large enough mass in stars exists that the IMF is sampled up to stellar masses large enough that their feedback begins to evaporate and/or disrupt the dense gas around them, and the SFR begins to decrease again.

The evolution of the SFR has been modeled by Zamora-Avilés, Vázquez-Semadeni & Colín (2012, hereafter ZA+12) and Zamora-Avilés & Vázquez-Semadeni (2014, hereafter ZV14; see also () and ()). The model, hereafter referred to as the ZV14 model, considers, similarly to other models of the SFR (Krumholz & McKee, 2005; Hennebelle & Chabrier, 2011; Padoan & Nordlund, 2011; Federrath & Klessen, 2012), that the turbulence generates a probability density distribution (PDF) of the density field, which in particular for nearly isothermal flows (such as molecular clouds) adopts a lognormal form (Vázquez-Semadeni, 1994; Passot & Vázquez-Semadeni, 1998), whose characteristic width depends on the turbulent Mach number (Passot & Vázquez-Semadeni, 1998; Federrath, Klessen & Schmidt, 2008). However, those other models assume that the turbulence opposes the gravitational contraction and is strong enough to prevent the global collapse of the clouds, thus being stationary in nature, and giving predictions, for example, for the star formation efficiency per free-fall time, ${\rm SFE_{\rm ff}}$, since the final efficiency depends on the lifetime of the clouds. Those stationary models then compute the instantaneous SFR in the cloud as the mass at high densities divided by a suitable free-fall time representative of those high densities. The various models differ in their choice of the “high densities” and their characteristic timescales.

Instead, the ZV14 model assumes that the turbulence is relatively weak (sonic Mach number of the order of a few, rather than of order 10, as is often assumed) and insufficient to support the cloud, so that the cloud can proceed to global and hierarchical gravitational contraction. Thus, this model can predict the temporal evolution of the SFR and the SFE, from which suitable temporal and ensemble averages can be computed that correctly match observations (ZV14). A brief summary of the ingredients of the model has been given by Vázquez-Semadeni et al. (2018).

The ZV14 model then assumes that, due to the global cloud contraction, the density PDF shifts to higher density, as the cloud’s mean density increases.555The model assumes that the turbulent motions proper remain at a roughly constant moderate level throughout the evolution, and that the strongly supersonic motions routinely observed in MCs are dominated by infalling motions that do not oppose the collapse, but rather result from it. This in turn implies that the mass at high densities, and consequently the SFR, initially increase over time. This is because, as the cloud contracts, a progressively larger fraction of its mass is at sufficiently high densities to undergo “instantaneous” collapse.

The model also considers a standard IMF to compute the mass in massive stars produced each timestep, given the total stellar mass formed up to that time. Next, from the total number of massive stars present at every time, the model calculates the photoionizing radiation produced, which begins to disrupt the cloud. So, eventually the cloud’s mass begins to decrease again.

It is important to note that the ZV14 model assumes that the initial conditions for the clouds are those of the cold atomic gas (e.g., Heiles & Troland, 2003) in the conditions under which it arises from nonlinear triggering of the thermal instability by converging flows in the warm medium that occupies most of the Galactic disk at Solar Galactocentric radii. Thus, contrary to the stationary models, the initial density and turbulent Mach numbers are fixed, rather than free parameters of the model. The only remaining parameter is the total mass of gas that becomes Jeans unstable and proceeds to collapse.

Figure 7 shows the evolution of clouds of masses $10^{3}$, $10^{4}$, $10^{5}$, and $10^{6}M_{\odot}$ according to the ZV14 model. It is important to note that, by construction, the ZV14 model matches the evolution of the SFR in the numerical simulations of Vázquez-Semadeni et al. (2010), implying that this evolution is representative of that observed in numerical simulations.

This kind of evolution has several important implications. First, the model predicts that young clouds will have low SFRs. This eliminates the supposed need for a supporting agent in order to keep low-mass star-forming clouds from having too large SFRs. The clouds are proceeding at free-fall, but, as discussed in Sec. 3.3, and indicated by eq. (11), this process takes time, and develops very slowly at its early stages (see also Burkert & Hartmann, 2013; Girichidis et al., 2014). As a consequence, the SFR is very low at early stages.

Indeed, ZA+12 calculated the evolutionary track of a 2000-$M_{\odot}$ cloud in the Kennicutt-Schmidt (KS) diagram of $\Sigma_{\rm SFR}$ vs. $\Sigma_{\rm gas}$ (Kennicutt, 1998). We show this track in Fig. 8, indicating the times before the final burst that destroys the cloud, together with observational data for various types of clouds, from low-mass star forming clouds (Evans et al., 2009; Lada et al., 2010, 2013) to high-mass regions like the OMC1 cloud and other mini-bursts (Louvet et al., 2014).

Figure 8 shows that, during its early stages, the cloud passes through the locus of low-mass star-forming clouds, and then, during its later stages, passes near the locus of high-mass star-forming clouds, such as the Orion OMC1 cloud, and finally, near the locus of mini-bursts, clearly indicating the increase in its star formation activity. It is important to remark that, as the SFR increases, the cloud shrinks, as indicated by the top-right panel of Fig. 7, so that the “cloud” is already a dense, compact, massive star-forming clump when it reaches SFRs characteristic of massive star-forming regions. Also, the clump at this time is expected to be embedded in a more massive cloud, as accretion is a fundamental part of the evolutionary process, and so clouds are expected to continue accreting from their surroundings during their evolution.

It is also important point to note that, according to the ZV14 model, all parameters of a cloud, such as its dense mass fraction, mean density, size, and SFR, evolve simultaneously. Therefore, pairs of parameters, such as the total instantaneous molecular mass and the dense mass fraction, can be used to completely determine the evolutionary state of a cloud of a given total dense gas mass. This property of the model has been used by Vázquez-Semadeni et al. (2018) to estimate the ages of the clouds in the sample studied by Lada et al. (2010), using the dense gas mass fraction as a proxy for evolutionary time. This study concluded, for example, that clouds of similar masses but very different instantaneous SFRs and SFEs, such as the California and the Orion B clouds, can be interpreted as being in very different evolutionary stages. For these two clouds, the California cloud is described by the model as being only a few megayears old, while the Orion B cloud is described as being between 15 and 20 Myr old.

Another point to remark is that ample observational evidence exists for the acceleration of star formation, either from the presence of tails of old-age stars in the age histograms in young, embedded clusters (Palla & Stahler, 1999, 2000; Da Rio et al., 2010), or by the more numerous and more centrally concentrated nature of the younger members of a cluster (e.g., Povich et al., 2016; Caldwell & Chang, 2018; Großschedl et al., 2019).

5.4.3 On the usage of the lognormal form of the PDF

It is important to emphasize that the ZV14 model assumes that the density PDF has a lognormal form, in spite of the fact that it is now well known that the PDF develops a power-law tail at high densities as various collapses proceed throughout the cloud (e.g., Kainulainen et al., 2009; Kritsuk et al., 2011; Ballesteros-Paredes et al., 2011b; Girichidis et al., 2014). A similar model using the growth of the power-law tail of the PDF has been presented by Burkhart (2018). However, it is not obvious that the power-law form is actually the most appropriate form to use for the purpose of the model. This is because, as suggested by Kritsuk et al. (2011), the power-law tail of the PDF appears to be related to the development of power-lar radial density profiles in collapsing regions. In this interpretation, the power-law part of the PDF represents the already-collapsing material in the cloud. Using the lognormal PDF attempts to represent the density fluctuations produced exlusively by the turbuelence. This is equivalent to assuming that the velocity field in a turbulent, self-gravitating cloud can be assumed to consist of a purely turbulent component plus an infall component, so that the infall velocities are associated with the power-law PDF, and the purely turbulent velocities are associated with the lognormal form, and that it is the latter that constitute the seeds for the subsequent collapse.

5.5 Combining fragmentation and the acceleration of star

formation. The assembly of clusters

Numerical simulations of the global contraction of clouds show that the first fragments to complete their collapse have lower masses than those that complete their collapse later. This is because the fluctuations that terminate their collapse much earlier than the bulk of the cloud are the ones that have much shorter free-fall times (i.e., much larger densities) than the mean values in the cloud. According to the lognormal form of the turbulent density PDF, the total mass in these high-density fluctuations is a small fraction of the total cloud mass, and therefore, these first collapses involve low-mass objects. This is illustrated in Fig. 9, which shows the evolution of the mass fraction at densities above $n_{\rm SF}=10^{3}\,n_{0}$, where $n_{\rm SF}$ is defined in Sec. 3.5, using the same approach as in Sec. 3.5. This mass fraction is seen to increase in time, so that initially it is very small, implying that the first star-forming sites within the cloud that appear must be of low mass. In turn, this suggests that the stars that form should also be of low mass. As time progresses, the mass in these sites increases, and thus they may be more massive themselves, as well as the stars they contain, assuming that each collapsing site forms a distribution of stellar masses bounded by the mass of the site itself (e.g., Oey, 2011).

The temporal increase of the mass involved in star formation within a cloud also implies a temporal increase of the SFR, as observed in numerical simulations (e.g., Heitsch & Hartmann, 2008; Vázquez-Semadeni et al., 2010; Hartmann et al., 2012; Colín et al., 2013; Körtgen et al., 2016), and predicted by the ZV14 model. Moreover, numerical simulations of whole GMCs (at the hundred-pc scale) in which the sink particles represent individual stars rather than stellar groups, also exhibit a temporal increase of the mean mass of the stellar particles. This is shown in Fig. 10, reproduced from Vázquez-Semadeni et al. (2017), where it is seen that the stellar-particle mass distribution evolves by developing a larger fraction of high-mass particles as time progresses, in tandem with the increase of the SFR. Finally, the formation of subsequently more massive stars while the cluster itself becomes more massive implies a correlation between the mass of the most massive star and that of the cluster itself, as is indeed observed (Weidner & Kroupa, 2006).

Once the SFR is high enough, the feedback from the massive stars that appear in the cloud begins to erode the cloud, and causes the SFR to decrease again, together with the high-mass fraction. This evolution of the SFR implies that the stellar age histograms of star-forming regions must peak at the age corresponding to the maximum of the SFR and to decrease toward both larger and smaller ages, very much in agreement with observed age histograms of embedded clusters (e.g., Palla & Stahler, 1999, 2000; Da Rio et al., 2010; Vázquez-Semadeni et al., 2017).

The spatial distribution of the stars in young clusters is also predicted to be affected by the GHC mechanism. Recall that filaments constitute the large-scale accretion flow (cf. Sec. 5.3), and that small-scale collapses are the first ones to terminate, but that they also accrete. Thus, the small-scale collapsing sites can continue to form stars as they fall into the large-scale potential well, feeding the latter with both stars and gas.

The stars that fall into the central hub inherit the velocity of the infalling gas. However, when both the stars and gas reach the hub, the gas is shocked (see Fig. 4 of Gómez & Vázquez-Semadeni, 2014), and loses kinetic energy. Therefore, stars formed later in the hub from the shocked gas, inherit this velocity dispersion, and thus have lower typical velocities than the stars that fell from the outside. Thus, the newer stars tend to have smaller orbits than the ones that fell in from the outside, generating an age gradient in the cluster, as shown in Figs. 8 and 9 of Vázquez-Semadeni et al. (2017). Recent young cluster studies suggest that such age gradients are indeed observed (e,g, Getman et al., 2014a, b).

Finally, not all the material in the filaments feeding the hub has time to reach it before it is destroyed by the feedback from the massive stars formed there. For example, the time for massive stars to appear after the formation of the first stars observed in the simulation studied by Vázquez-Semadeni et al. (2017) is $\sim 5$ Myr. After that, an Hii region begins to expand, with a characteristic speed $\sim 10$-$20\,{\rm~{}km~{}s}^{-1}$. On the other hand, material in the filaments, at least for the range of masses of the structures formed in the simulation, has typical infall speeds $\sim$ a few${\rm~{}km~{}s}^{-1}$. Therefore, material in the filaments located at distances greater than $\sim 5$-10 pc from the main hub is likely to be exposed to the ionizing radiation in the Hii region before it can fall into the hub. This radiation evaporates first the filament, interrupting the supply of gas to the cores within it (or “starving” them), leaving for a while a trail of dense cores embedded in the warm ionized gas, similar to Bok globules. Later, the cores are consumed by their internal star formation and/or by evaporation by the ambient radiation. This is illustrated in Fig. 11, which shows a filamentary cluster in the simulation LAF1 of Colín et al. (2013) at 1.6 and 8.6 Myr after the formation of the first star. The Hii region begins to expand some $\sim 5$ Myr after the formation of the first star. The filament in the lower part of the images is seen to be evaporated, leaving behind a chain of isolated dense cores (“globules”). The cores are finally consumed too, roughly 1 Myr later. A similar behavior is observed in the simulations including full radiative transfer of Zamora-Avilés et al. (2019).

If the distant cores along the filaments have already formed some stars, then, when they are consumed/evaporated, the stellar groups formed there exhibit a hierarchical and fractal structure similar to that established originally in the gas, as shown in Fig. 12 of Vázquez-Semadeni et al. (2017).

6 Discussion

In the previous sections we have outlined the mechanism of GHC and fragmentation, extending the scenario originally proposed by Hoyle (1953) to the case of clouds containing nonlinear density fluctuations produced by (moderately) supersonic turbulence, as motivated by the observed evolution of numerical simulations of the formation and self-consistent turbulence generation and subsequent GHC. In what follows, we discuss some important insights and implications of the scenario, as well as some caveats of our approach.

6.1 The IMF: Minimum fragment mass and the formation of brown dwarfs

The minimum mass that can be reached through the GHC mechanism is an important question that needs to be addressed in order to understand whether brown dwarfs can be produced by this mechanism. Already in his original paper, Hoyle (1953) gave an estimate for this minimum mass, assuming that it would be determined by the moment at which the collapsing gas becomes dense enough to become optically thick, and thus suffer a transition from a nearly isothermal to a nearly adiabatic regime. From this reasoning, he concluded that the minimum mass should be of order $\sim 0.3$–1.5 $M_{\odot}$. Further refinements on the cooling mechanisms yielded significantly lower masses, down to $\sim 0.01\,M_{\odot}$ (Rees, 1976, and references therein), while refinements on the assumed geometry, such that the fragmentation occurs in a shock-compressed layer and that accretion causes a fragment to grow while still condensing out, lead to more modern estimates of the minimum mass $\sim 3\times 10^{-3}\,M_{\odot}$ (Boyd & Whitworth, 2005).

Numerical simulations including this so-called opacity-limited fragmentation also show an efficient formation of brown dwarfs (Bate et al., 2002), and in fact the problem then becomes to avoid an excessive formation of these objects by inclusion of radiative transfer to adequately model the gas heating from the accretion shock onto the hydrostatic core, in order to prevent excessive fragmentation (e.g., Krumholz et al., 2007; Bate, 2009b; Offner et al., 2009). Although these simulations refer to fragmentation in cores within initially turbulence-supported clouds according to the “gravoturbulent” scenario, the last stages of fragmentation must be very similar as in our scenario, so their conclusions also apply to the GHC scenario. Therefore, it is safe to conclude that the GHC fragmentation scenario can readily produce objects down to the brown dwarf regime. See Sec. 6.7 for a comparison between GHC and the gravoturbulent scenario.

6.2 The IMF: the upper-end as a consequence of chaotic GHC

The GHC scenario also provides a natural framework for the development of the Salpeter (1955) slope of the stellar initial mass function (IMF). The accretion onto stars, or protostellar objects, constitutes the last step of the mass (accretion) cascade that must start at the scale of the contracting cloud. Numerical simulations using a variety of schemes and setups have successfully reproduced the high mass end of the IMF, with a slope $\sim-1$, close to the $-1.3$ canonical value of Salpeter (e.g., Bate, 2009a; Bonnell et al., 2011). The usual approach in those simulations has been to create sink particles—particles that continue interacting gravitationally, but not hydrodynamically, with their environment— in places where the conditions are believed to be adequate for star formation, such as large densities, gravitational binding against all possible sources of support (magnetic fields, turbulence, thermal pressure), and local inward motions. Once created, the sink particles are allowed to continue accreting, increasing their mass with time, according to the local conditions (e.g., Bate et al., 1995; Jappsen et al., 2005; Federrath et al., 2010). In principle, when the simulation is properly resolved, and the relevant physics is included, the distribution of sink-particle masses should reproduce the IMF, since the accretion onto the sink particles should mimic the actual physical process occurring in stars. In the rest of this section we discuss how the GHC scenario allows the accretion processes that may lead to the Salpeter slope. We do not discuss the turnover of the IMF, which may be more related to local stellar physics than to the global gravity of the medium (e.g., Krumholz et al., 2007; Bate, 2009b).

6.2.1 The kinetic or Boltzmann approach. The last stage of

accretion

While the thermal and kinetic physics of the numerical simulations that have obtained a sink mass function with slopes close to the Salpeter value varies from one work to another, all of them have one common ingredient: gravity. In particular, Zinnecker (1982) proposed an approach to the Salpeter slope based on gravitational accretion assuming that the process of accretion of mass onto the protostellar objects is of a Bondi-Hoyle-Littleton (BHL; Hoyle & Lyttleton, 1939; Bondi, 1952) type. In this scheme, a star of mass $M$ traveling with a relative velocity $v$ through a homogeneous medium of density $\rho_{0}$ presents an effective cross section determined by the star’s gravitational potential well and its velocity relative to the medium. The resulting mass accretion rate is

[TABLE]

with $c_{s}$ the sound speed of the medium. Zinnecker (1982) showed that, for a given population of protostellar objects and $\alpha$ constant, the mass distribution of that population will evolve to an asymptotic power-law slope of

[TABLE]

where $\Gamma$ is the exponent in the logarithmic mass distribution of the stars, i.e.,

[TABLE]

and $dN/dM\propto M^{\Gamma-1}=M^{-2}$.

Although this approach gives a slope $\Gamma=-1$ close to the Salpeter value, it nevertheless has some problems:

The environment where the star is being formed is far from homogeneous. 2. 2.

It neglects the self-gravity of the medium. 3. 3.

The solution to eq. (26) diverges in a finite time.

Indeed, numerical simulations (e.g., Klessen & Burkert, 2000; Bonnell et al., 2006; Maschberger et al., 2014; Ballesteros-Paredes et al., 2015) exhibit departures from the pure BHL accretion:

The mass accretion histories of the sink particles have time-variable slopes, starting very steep, and then decreasing to being almost flat in most cases (see, e.g., Fig. 13 of Ballesteros-Paredes et al., 2015). 2. 2.

The ensemble of stars in simulations does not exhibit a clear $\dot{M}\propto M^{2}$ accretion law, but instead exhibits a scattered distribution in the $\dot{M}$–$M$ plane, the scatter being larger at lower $M$ (see, e.g., Fig. 2 of Ballesteros-Paredes et al., 2015). The lower boundary of the distribution has a slope $\sim 2$, but at the centermost parts of the distribution, it has a flatter slope. 3. 3.

The individual mass accretion rates exhibit large oscillations during the accretion process, likely due to the fluctuations in the environmental conditions (Maschberger et al., 2014; Kuznetsova et al., 2018b).

Even though these points suggest that BHL is not a complete model for the mass accretion onto the stars, it comes reasonably close. Indeed, as shown by Ballesteros-Paredes et al. (2015), the dependence $\dot{M}\propto M^{2}$ is a good approximation if objects are collected into groups with similar values of $\alpha$; i.e., groups of objects forming in regions where the density and velocity dispersion are similar. This is shown in Fig. 12, where the mass accretion rate divided by the local $\alpha$ is plotted against the mass of the sinks, in simulations of GHC, after one free-fall time of the initial density. Each range in $\alpha$ is denoted with a different symbol and color. The dashed line has a slope of 2, which is seen to be close to the slope defined by objects depicted with the same color/symbol.

In view of this, Ballesteros-Paredes et al. (2015) argue that, since the environment of the sinks is temporally and spatially variable, one may write the accretion rate in a more general form as

[TABLE]

and therefore, that Zinnecker (1982)’s approach is applicable even if $\alpha$ is not constant, as long as it does not depend on $M$. This generalization of the formalism to describe environmental variability can account for the mass accretion rates observed in simulations, in comparison to those predicted by the BHL approach.

6.2.2 The core mass function and the IMF: the

next-to-last accretion stage

It is well known that the mass function of the dense cores that form stars (core mass function, or CMF), bears a striking similarity to the stellar IMF (Simpson et al., 2008; André et al., 2010), which has often been interpreted as evidence that the final stellar masses are determined already at the dense core stage (e.g., André et al., 2014). The connection between the slope of the core mass function (CMF) and the IMF becomes natural in the global-collapse scenario, since one of its essential features is that all scales accrete from their parent structures. Indeed, as shown by Kuznetsova et al. (2017), $n-$body simulations of cluster formation by gravitational contraction without the effects of gas physics, develop both the quadratic dependence of the mass accretion rate on mass, $\dot{M}\propto M^{2}$ and the mass function of clusters with slope $\Gamma=-1$. This is in agreement with the well-known result that the cluster mass function also exhibits a slope of $-1$ (see, e.g., Oey, 2011, and references therein). More relevant for our interests here, is the fact that a similar result is found by Kuznetsova et al. (2018b) when the gas physics is included, suggesting that the effect is purely gravitational. This suggests that the BHL approach is applicable, as a first approximation, to the accretion onto the dense cores where the stars form, as well as to the stars themselves.

Thus, in the GHC scenario, the resemblance between the CMF and the IMF may be the result of a similar accretion mechanism operating both at the core scale and at the protostellar object scale. This is in contrast to “core collapse” models (e.g., Padoan & Nordlund, 2002; McKee & Tan, 2003; Hennebelle & Chabrier, 2008; Hopkins, 2012) where the star simply “inherits” a fraction of the the mass of its parent core. While compact pre-stellar cores may be gravity-dominated, they still do fragment into several smaller units, frequently simply termed “fragments” (e.g., Palau et al., 2014, 2015; K. Lee et al., 2015; Beuther et al., 2018). But both the cores and the stars will tend to have a Salpeter-like mass distribution as long as their accretion process is dominated by the local gravity.

6.2.3 The IMF in strongly turbulent environments

In “core collapse” or “gravoturbulent” models of star formation and the IMF (e.g., Padoan & Nordlund, 2002; McKee & Tan, 2003; Hennebelle & Chabrier, 2008; Hopkins, 2012), the CMF is assumed to be the result of pure turbulent fragmentation, and the cores then simply contract to form the stars, while the cloud continues to be globally supported by turbulence. However, the results of several numerical simulations (Clark et al., 2008; Bertelli Motta et al., 2016) suggest instead that, when turbulence dominates (i.e., for relatively large values of the virial parameter, $\alpha_{\rm vir}\ga 3$), the number of low-mass stars decreases, the mass distribution becomes more top-heavy, and thus, the slope of the mass function of sinks departs from the Salpeter slope. The reason for this behavior is that in the presence of strong turbulence, only high mass pre-sink entities have enough gravitational energy to overcome the kinetic energy and proceed to collapse. We conclude that, rather than turbulence, what sets the masses of cores and stars is a BHL-type accretion mechanism. Stars may have a slightly steeper slope because the higher masses are limited by the mass of the parent clump where they form (Oey, 2011). On the other hand, uncertainties in the measured IMFs are often large enough that it may be impossible to distinguish between a slope of $-1$ and a slope of $-1.3$. In fact, slopes of $-1$ (or $-2$ in $dN/dM$) are often reported (e.g., De Marchi et al., 2010).

6.3 Simultaneous evolution of cloud physical properties and the

star formation activity

An essential feature of the GHC scenario is that it is evolutionary at the cloud level (see further discussion in Sec. 6.7). The evolution occurs at the level of all the cloud properties: their masses, densities, dense mass fractions, energy balance, and star formation activity. In particular, in Sec. 5.4.2 we discussed the increase of the SFR in a cloud due to the increase of its mean density and, as a consequence, its dense gas fraction. This result offers a straightforward explanation for the observed correlation between the fraction of clumps associated with massive stars and their peak column density (e.g., Urquhart et al., 2018) or between the surface density of dense gas and of the star formation rate (e.g., Gao & Solomon, 2004; Lada et al., 2010), as shown by Vázquez-Semadeni et al. (2018, see also Secs. 5.4.2 and 6.6.2, and Camacho et al. 2019, in prep.).

6.4 Evolution of the energy balance of clumps and cores

As discussed by Ballesteros-Paredes et al. (2011a) and in Sec. 3.2, gravitationally contracting objects are expected to attain contraction velocities of the order of the free-fall, or gravitational, speed, eq. (4) after roughly one free-fall time (cf. Sec. 3.2). Because this speed differs by only a factor of $\sqrt{2}$ from the “virial” speed

[TABLE]

the observation that clouds and their substructures tend to appear “virialized” can be understood if these structures are contracting at approximately the gravitational speed.

When the velocity is originated by the gravitational contraction of the objects, it is of the order of the gravitational speed, eq. (4). In this case, a scaling relation between the Larson ratio ${\cal L}$, (cf. Sec. 2.2) and the column density of the form

[TABLE]

is expected, and indeed is observed in general (e.g., Keto & Myers, 1986; Heyer et al., 2009; Ballesteros-Paredes et al., 2011a; Field et al., 2011; Leroy et al., 2015). We refer to this condition generically as equipartition between the kinetic and self-gravitating energies.

6.4.1 The approach to equipartition and the inverse

$\alpha_{\rm vir}$-mass relation

Although in general most MCs and their substructures appear to be near the equipartition relation, eq. (31), significant deviations are routinely observed around this scaling. In particular, an inverse $\alpha_{\rm vir}$-mass correlation is often observed in surveys ranging from GMCs (e.g., Keto & Myers, 1986; Leroy et al., 2015; Traficante et al., 2018a) to dense molecular cores (e.g., Kauffmann et al., 2013; Ohashi et al., 2016; Sanhueza et al., 2017; Traficante et al., 2018b). Because of the variety of objects considered in this section, we will generically refer to them simply as “the objects”. This inverse correlation can be understood within the context of the GHC scenario, in combination with a common selection effect inherently introduced by the survey observational procedure.

A fundamental point to understand the $\alpha_{\rm vir}$-mass correlation is that the turbulent (or “inertial”) motions, characterized by a velocity dispersion $\sigma_{\rm turb}$, are independent of the self-gravity-driven motions in the objects, which in turn are characterized by the gravitational velocity $v_{\rm g}$ (see Sec. 3.2, and eq. [5] therein). This consideration differs essentially from the standard assumption that the turbulence somehow tend to be “virialized”, an assumption that we have questioned in Sec. 2.2. Indeed, throughout this paper we assume that the turbulent motions do not automatically adjust to the local gravitational potential to “virialize”, but rather constitute just a background on top of which the self-gravitational motions develop, and create a tendency toward equipartition.

This coexistence of the turbulent and gravitational motions is incorporated into the virial parameter in eq. (9). Combining this equation with eq. (10), we can obtain an equation for the prestellar evolution of the virial parameter of a gravitationally contracting object of constant666Recall the derivation of eqs. (8) and (9) in Sec. 3.2 was made assuming a constant mass of the object. This treatment neglects the accretion onto the object, and has the consequence that the predicted trajectories of objects in the ${\cal L}$-$\Sigma$ and $\alpha_{\rm vir}$-mass diagrams differs somewhat from those oberved in simulations (Ballesteros-Paredes et al., 2018, see also Camacho et al. 2019, in prep.). We discuss this some more in Sec. 6.10. mass $M_{i}$, where the subindex $i$ denotes the $i$-th object of a survey:

[TABLE]

where $t$ is the time since the object began contracting, and we have used the facts that $\tau_{i}=t/t_{{\rm ff},i}$ and $t_{{\rm ff},i}=\sqrt{3\pi/32G\rho_{0,i}}$, where $\rho_{0,i}$ is the density of the object when it began contracting. Note that the first term in the right-hand side (RHS) of eq. (32) can be either smaller or larger than unity, while the remainder of the RHS varies between 0 and 1 for $0\leq t\leq t_{\rm ff}$.

Note also that, due to the standard definition of the virial parameter, the denominator contains the final free-fall speed $v_{\rm ff}$, given by eq. (4), rather than the actual (smaller) gravitational speed $v_{\rm g}$, given by eq. (5). The first term within the curly brackets then reflects the contribution of the turbulent (or “inertial”, as opposed to driven by self gravity) velocity to $\alpha_{\rm vir}$, normalized by the final free-fall speed, while the rest of the right hand side reflects the contribution of the true self-gravity-driven speed to $\alpha_{\rm vir}$, and shows that a given time $t$ represents a larger fraction of the free-fall time for denser objects, because their free-fall time is shorter. Therefore, denser objects will be closer to the final, free-fall speed than lower-density objects of the same age.

For an individual object, eq. (32) shows that, as it contracts gravitationally, it should approach a value $\alpha_{\rm vir}\sim 2$. This is because, as the object contracts at constant $M$ (cf. Sec. 3.2), its $v_{\rm ff}$ increases, and so the ratio $\sigma_{\rm turb}^{2}/v_{\rm ff}^{2}$ decreases (recall we assume that the turbulent speed is independent of the gravitational potential). Also, as $t\rightarrow t_{\rm ff}$, the term within square brackets in the first equality of eq. (32) approaches zero, and so $\alpha_{\rm vir}\rightarrow 2$. However, it can approach it from above of from below, depending on the inital value of the ratio $\sigma_{\rm turb}^{2}/v_{\rm ff}^{2}$, as described in Sec. 3.2 and Fig. 2: If initially $\sigma_{\rm turb}^{2}>v_{\rm ff}^{2}$, then initially $\alpha_{\rm vir}>2$, and $\alpha_{\rm vir}\rightarrow 2$ from above. Conversely, if initially $\sigma_{\rm turb}^{2}\ll v_{\rm ff}^{2}$, then $\alpha_{\rm vir}\rightarrow 2$ from below.

In addition, eq. (32) can also be used to describe a sample of roughly coeval777I.e., objects that began contracting within a time span significantly shorter than their free-fall times. objects at a single moment in time, such as the dense cores within a clump, or fragments within a core, whose onset of collapse may be partially synchronized by the global reduction of their parent object’s Jeans mass, as decribed in Sec. 3.3. In this case the subindex “$i$” denotes the $i$-th object in the sample. In this case, however, in order to relate the objects’ densities to their masses, the selection effects introduced by the observational procedure must be considered. It is often the case that observational samples performed with a given instrument tend to define the sample members such that their mass scales as some power of their size, $M\propto R^{p}$, where typically $1\la p\la 3$ (e.g., Kauffmann et al., 2010, 2013; Simpson et al., 2011; Ragan et al., 2013; Wienen et al., 2015; Pfalzner et al., 2016; Traficante et al., 2018b).888Slopes $p\ga 2$ can be understood as the result of the observations requiring a certain minimum column density in order to be detectable above the noise, which effectively defines the objects through column density thresholds. In this case, Ballesteros-Paredes et al. (2012) showed that the mean column density of the sample members is close to the detectability threshold, because the typically steep slope of the column density PDFs of molecular clouds and their substructures implies that most of the material is at the lowest densities. On the other hand, slopes $1\la p\la 2$ typically arise when the object samples are created by using multiple thresholds, as in the case of multi-tracer observations (see Sec. 4.1. of Ballesteros-Paredes et al., 2019)

Taking into account this selection effect, we can consider the contributions from the inertial and gravitational terms in eq. (32). Concerning the inertial term, $\sigma_{\rm turb}^{2}/v_{{\rm ff},i}^{2}$, we note that $v_{{\rm ff},i}$, given by eq. (4) for object $i$, is larger for larger $M_{i}$ provided that $p>1$. But in addition, we must consider the size dependence of the turbulent velocity dispersion, given by exponent $\eta$ in eq. (7). Plausible values of $\eta$ are zero if the turbulent background is the same for all members of the sample, $1/3$ if the turbulence is subsonic and thus follows a Kolmogorov scaling, and $1/2$ if the turbulence is supersonic and follows a Burgers scaling. Therefore, the ratio $\sigma_{\rm turb}^{2}/v_{{\rm ff},i}^{2}$ is smaller for larger $M_{i}$ when $p>1$, 1.67 or 2 for $\eta=0$, 1/3 or 1/2, respectively.

Regarding the term within square brackets in eq. (32), we can find its dependence on $M_{i}$ as follows. If we write in general

[TABLE]

for the $i$-th member of the sample, then $M_{i}\propto\rho_{i}^{p/(p-3)}$, and so $M_{i}$ varies inversely with $\rho_{i}$ for $0<p<3$. Note that, here, $\rho_{i}$ is the density if the $i$-th object at time $t$, while eq. (32) involves $\rho_{0,i}$, the density of object $i$ when it began contracting. However, it is clear that, if $\rho_{i}(t)>\rho_{j}(t)$ for two objects $i$ and $j$ of the coeval sample, then $\rho_{0,i}>\rho_{0,j}$ as well, because the denser object increases its density faster than the less dense one. Thus, more massive objects within a sample that satisfies $M_{i}\propto R_{i}^{p}$ with $0<p<3$ have lower values of $\rho_{0}$. This implies that the term within the square brackets in the second equality of eq. (32) is smaller for the more massive objects of a roughly coeval sample, if the sample satisfies eq. (33) with $0<p<3$.

These results for a sample of objects can be applied to specific types of objects, as we now discuss.

Application to diffuse clouds and clumps

A clear tendency towards exhibiting values of the virial parameter $\alpha_{\rm vir}$ or the Larson ratio ${\cal L}$ significantly larger than their respective equipartition values is often reported in surveys of relatively low-column density clouds and clumps (e.g., Keto & Myers, 1986; Leroy et al., 2015; Traficante et al., 2018a), implying an excess of kinetic energy over the energy from the self-gravity of the cloud.

This kinetic energy excess can be interpreted in terms of eq. (32) mostly as relating to the first term, since for diffuse objects, which are also generally rather large (of sizes of parsecs and larger), the gravitational velocity is expected to start out smaller than the inertial velocity, and so the term $\sigma_{\rm turb}^{2}/v_{{\rm ff},i}^{2}$ dominates $\alpha_{\rm vir}$. Note, however, that this term being significantly larger than unity does not necessarily imply turbulent support or disruption. Any external compression that triggers the formation of a self-gravity-dominated object must start with the inertial compression velocity being larger than the gravitational velocity, or otherwise the object would have started being self-gravitating already. Thus, the formation of self-gravitating objects by converging flows or “collect and collapse” processes must start with an apparently strongly supervirial object.

This mechanism has been investigated by Camacho et al. (2016) in a survey of clouds, clumps and cores (generically, “clumps”) in numerical simulations of cloud formation and evolution by converging motions in the WNM. Those authors found that, in objects exhibiting such an excess in kinetic over self-gravitational energy, in roughly half of the cases the excess kinetic energy was in net compressive motions, while in the rest, the clumps were truly in the process of dispersal. That is, contrary to suggestions that these objects are in equilibrium, confined by large external pressures (Keto & Myers, 1986; Field et al., 2011), Camacho et al. (2016, see also Ballesteros-Paredes et al. 1999a) showed that the objects are not in equilibrium. Instead, they are either in the process of growth, by external compressions whose origin is other than their self-gravity, or in the process of dispersal. There is no reason to assume that such structures should in general be in equilibrium, and therefore, there is no need for a large confining thermal pressure.

These initially super-virial structures are then expected to evolve according to the process described by the dash-triple dot lines in Fig. 2 (if they are defined in a Lagrangian way, with a fixed mass, as in Sec. 3.2): an external compression (due to turbulence in the ISM, or a large-scale instability in the Galactic disk) induces the formation of the cloud, which, initially, is not dominated by its self-gravity. Thus, the cloud appears super-virial in the ${\cal L}$-$\Sigma$ or $\alpha_{\rm vir}$-$M$ diagrams. However, as its density increases due to the compression, and perhaps to thermal instability, the self gravity of the cloud increases, causing its gravitational speed $v_{\rm g}$ (eq. [5]) and Larson ratio to also increase, and its $\alpha_{\rm vir}$ to decrease. Thus, the energy in the self-gravitating velocity eventually becomes larger than the inertial compressive energy. At this point, the cloud becomes dominated by its self-gravity. Nevertheless, the cloud may then appear somewhat subvirial, if the inertial compressive motion dissipates rapidly and $v_{\rm g}$ has not yet become close enough to the free-fall value $v_{\rm ff}$ (eq. [4]) (the lower dash-triple dot curve in Fig. 2). Numerical simulations confirm the general shape of the evolutionary tracks of clumps in the ${\cal L}$-$\Sigma$ diagram (see the solid curves in the panels labeled “Mach 16” and “Mach 8” in Fig. 4 of Ballesteros-Paredes et al., 2018). Camacho et al. (2019, in prep.) show, in addition, the simultaneous evolution of the clumps’ energy balance and of their SFR.

Note that in no case it is necessary to invoke a large external confining thermal pressure for these apparently super-virial structures. If anything, the inflow, either inertial or from self-gravity, provides a ram pressure for the inner parts of the structures, although not really “confining” them, since the structures are not in equilibrium, but rather contracting. In the case of clumps that do not manage to become self-gravitating and collapse, then no confining pressure is necessary, either. They will simply disperse after the transient compression that formed them subsides.

Application to subvirial dense cores

On the other hand, dense cores are often observed to have values of ${\cal L}$ and $\alpha_{\rm vir}$ smaller than the gravitational value (e.g., Kauffmann et al., 2013; Ohashi et al., 2016; Sanhueza et al., 2017; Traficante et al., 2018b). The energy balance of these objects can also be understood in terms of the onset of their own collapse at a specific moment in time and with a specified initial radius, as described by eq. (8), if in this case the initial inertial velocity is smaller than the object’s own free-fall speed. In this case, the evolution of the object is generically described by the trajectories depicted by the solid or dashed curves in Fig. 2. It is important to point out that, for dense cores that are located within larger, less dense clumps which are already engaged in gravitational contraction, the onset of the cores’ own gravitational contraction is triggered by the global temporal decrease of the average Jeans mass in the clump, rather than by a strong transient, local reduction of the Jeans mass triggered by the inertial compression (Clark & Bonnell, 2005). But because the contraction begins with a finite radius and at a specific time, the infall speed is also smaller than the free-fall speed, and the core appears sub-virial, even though the core is nevertheless proceeding to collapse freely. Thus, its infall speed will asymptotically approach the free-fall speed as the core evolves. Again, numerical simulations confirm this generic form of the evolutionary track of cores for which the initial inertial motions are less than their own free-fall speed (see panels “Mach 4” and “Mach 0” in Fig. 4 of Ballesteros-Paredes et al., 2018).

It is important to note that, in this case, it is not necessary to invoke a strong magnetic field to support the cores. Their evolution will self-consistently make them appear sub-virial during the early stages of their collapse. Moreover, since the early stages of the collapse occur rather slowly, the cores will tend to spend more time in this sub-virial state than in the full-equipartition stage, making this subvirial state a common feature of these cores.

Note that the mechanism described in this section is independent of the recent suggestion by Traficante et al. (2018c) that the observed subvirial nature of massive cores may be due to an observational bias which causes the kinetic energy to be systematically underestimated in observations of such cores, due to the different tracers used to estimate the mass and the velocity dispersion. Although this effect may certainly be occurring in several massive core surveys, it is unlikely that it is the only explanation, since not all surveys employ the same tracers. Moreover, this effect cannot explain the observation that dense cores start out subvirial in numerical simulations (Ballesteros-Paredes et al., 2018, see also Camacho et al. 2019, in prep.), since in those numerical studies the energy determinations have been made directly from the density and velocity fields in the cores, with no intermediate synthetic observation step that could introduce the bias. Thus, we consider that the bias discussed by Traficante et al. (2018c) can be at play in some observational studies, but that the subviriality is a real physical property of most dense cores.

6.5 Outside-in vs. inside-out collapse in prestellar cores

Prestellar cores are often observed to 1) have Bonnor-Ebert-like (BE) column density profiles (e.g., Teixeira et al., 2005); 2) have molecular line profiles with a central self-absorption dip and a blue-peak excess that seem to imply subsonic infall speeds (e.g., Zhou, 1992; Evans, 1999; Lee et al., 2001; Campbell et al., 2016), and 3) to exhibit “extended inwards motions”, by which it is meant that they seem to extend beyond the expected radial location of the “rarefaction front” according to Shu’s inside-out collapse model (Lee et al., 2001). The first two of these properties are often interpreted as evidence of quasi-static contraction in the cores, as proposed by Shu (1977, herafter S77), although it has been difficult to reconcile them with the third property (e.g., Lee et al., 2001; Bergin & Tafalla, 2007). Moreover, low-mass cores often appear to be gravitationally unbound, thus requiring an external confning pressure to keep them from dispersing (e.g., Lada et al., 2008). In this section we discuss how these properties can be understood in the framework of the GHC scenario.

6.5.1 Dynamic vs. quasistatic prestellar contraction

The non-homologous and outside-in nature of the prestellar core contraction discussed in Sec. 4 allows a reinterpretation of the observational data. First and foremost, the fact that dynamical contraction starts much earlier than the formation of the first singularity (i.e., since the prestellar stage) is opposite to the assumption by S77 that the prestellar stage occurs quasi-statically, and that dynamical collapse begins at the time of the formation of the protostar (the singularity), leading to the well-known scenario of an inside-out collapse.

S77 proposed that the prestellar stages should contract quasi-statically rather than dynamically on the basis of two main arguments: First, that the conditions necessary for establishing a Larson-Penston (Larson, 1969; Penston, 1969, hereafter, LP flow) dynamical flow are ad hoc and difficult to establish in reality. Second, that the establishment of an $r^{-2}$ density profile represents the tendency of an isothermal, self-gravitating gas to approach detailed mechanical balance, and that this can be accomplished as long as different parts of the cloud can communicate acoustically with each other, which requires a subsonic flow. Measurement of apparently subsonic inflow speeds from infall line profiles (e.g., Evans, 1999, and references therein) have provided support to this view.

However, subsequent work has demonstrated the inapplicability of these arguments. Regarding the first of these, numerical simulations and analytical studies alike have long suggested that a wide variety of initial conditions generate flow that asymptotically approaches LP flow (e.g., Hunter, 1977; Whitworth & Summers, 1985; Foster & Chevalier, 1993; Naranjo-Romero et al., 2015). So, rather than being and artificial and ad-hoc condition, the LP flow appears to be an attractor for the collapse flow.

Concerning S77’s second argument, various pieces of evidence suggest that, rather than representing detailed mechanical balance, an $r^{-2}$ profile may simply be the manifestation of nearly pressureless collapse. Indeed, this profile arises in the outer regions of numerical simulations of spherical collapse. Far from the center, the internal mass is much larger than the mean Jeans mass (Naranjo-Romero et al., 2015). Moreover, it has recently been shown that this density profile can arise simply as a consequence of pressure-free collapse, under the condition that the infall speed is generated by self-gravity and that the mass inflow rate is independent of radius (Li, 2018), providing an alternative to the suggestion by S77 that it must arise from detailed mechanical balance, which in turn would result from acoustic (subsonic) coupling throughout the core. Thus, the $r^{-2}$ density profile in the outer parts of prestellar cores can be arrived at through dynamical collapse during the prestellar stage, rather than requiring a quasistatic process, as suggested by S77.

6.5.2 Dynamical contraction with Bonnor-Ebert-like density

profiles in prestellar cores

The dynamical contraction during the prestellar stage is also known to produce density structures that resemble a Bonnor-Ebert (BE) profile, with a flat central part and an $r^{-2}$ scaling in the outer parts (the envelope), as denoted in eq. (24). However, contrary to true BE spheres, which are hydrostatic equilibrium configurations, contracting prestellar cores have a non-zero infall speed at all radii, except at the core center, as indicated by eq. (24). Moreover, the decrease of the infall speed towards the core’s center (i.e., the outside-in nature of the profile) implies that the densest parts of the core do not have very large infall speeds during the prestellar stage, while the largest speeds occur at radii where the density is already decreasing, and therefore those speeds are downweighted in line profiles, giving the appearance that the cores have smaller infall speeds than they actually do (Loughnane et al., 2018). This may help reconcile the supersonic nature of the actual infall speed with the apparently subsonic values often derived from blue-skewed molecular line profiles (Zhou, 1992; Evans, 1999; Lee et al., 2001; Campbell et al., 2016), as well as explain the observed BE-like column density profiles in spite of the configurations not being in equilibrium.

In addition, if the cores have been undergoing dynamical collapse ever since the time when the local density fluctuation became unstable, the radial extent of the contraction motions must be much larger than the position of the rarefaction front in S77’s model, which only starts to propagate at the time of the formation of the protostar. That is, the local collapse motions have been propagating outwards for the whole prestellar stage of the contraction, while in S77’s model they only begin to propagate at the time of the formation of the protostar. This may explain the observation of “extended inwards motions”, at radii larger than the expected position of the rarefaction front in the inside-out model (Lee et al., 2001).

Finally, in this scenario, the boundary of cores that started as finite-extent turbulent density fluctuations may be defined in an observationally-motivated way, as the radius at which the density fluctuation merges into the background.999 Observationally, the background is often defined as the radius at which the signal-to-noise has become too low, or where the gradient of column density suffers an abrupt change (see, e.g., André et al., 2014). However, because the fluctuation grows by developing an $r^{-2}$ density profile and increasing its central density, the boundary of the core moves outwards as the core evolves. Thus, cores defined in this way grow both in size and mass. Naranjo-Romero et al. (2015) showed that, when the boundary of a collapsing core is defined in this way, the core’s evolution tracks the locus of observed cores in The Pipe and the Orion clouds in a diagram of $M_{\rm c}/M_{\rm BE}$ vs. $M_{\rm c}$, where $M_{\rm c}$ is the mass of the core, and $M_{\rm BE}$ is the BE mass for the average density and temperature of the core. This diagram was first investigated by Lada et al. (2008) for the cores in The Pipe. Those authors found that both stable ($M_{\rm c}/M_{\rm BE}<1$) and unstable ($M_{\rm c}/M_{\rm BE}>1$) cores occupied one common locus in this diagram. However, they concluded that the stable cores are gravitationally unbound and have to be confined by external pressure, while the unstable cores would be out of equilibrium, and collapsing. This interpretation, however, leaves the questions open as to why would both quiescent and dynamic cores occupy the same curve in this diagram. In fact, Lada et al. (2008) mention that the pressure on the cores is most likely due to the weight of the surrounding MC. This is consistent with the cores just being the densest, inner part of a large-scale collapsing object, where the infall speed is low, as dictated by eqs. (24), and the pressure being provided by the ram pressure of the external infalling material, as suggested by Heitsch et al. (2009).

Within this context, the suggestion by Naranjo-Romero et al. (2015) is that all cores in Lada et al. (2008)’s sample are collapsing, and that the “stable” ones only appear so because they have been observationally truncated at radii that are too short compared with the extension of the infall motions. Therefore, it would be highly desirable for observations of dense cores, to obtain data, when possible, up to distances as large as the extension of the inwards motions (e.g., Lee et al., 2001). Also, Naranjo-Romero et al. (2015) suggested that the observed BE-like shape of observed prestellar cores (e.g., Teixeira et al., 2005) is a manifestation of their being in a prestellar stage, and not indicative of hydrostatic equilibrium.

Finally, the protostellar stage (i.e., after the appearance of a singularity, or protostar) has been investigated analytically by Murray & Chang (2015) and numerically by Murray et al. (2017), who also conclude that the cores never go through a hydrostatic stage.

6.5.3 The decrease of the linewidth towards the central parts

of cores

In Sec. 4.2 we pointed out that the radial velocity profile in in the inner part of prestellar cores is of the form $v(r)\propto-r$ (cf. eq. 24) (e.g., Whitworth & Summers, 1985), implying that the infall speed decreases towards the center. This may offer a different interpretation for the observation that the nonthermal component decreases towards the core centers (e.g., Goodman et al., 1998; Pineda et al., 2010; Chen et al., 2018; Sokolov et al., 2018). This decrease is sometimes inferred from radial scans of the line profile moving away from the core’s center, in which a progressive blurring of groups of hyperfine lines is observed towards the outer regions of the core (e.g., Pineda et al., 2010). This blurring is assumed to be due to an increase in the turbulent component of the velocity dispersion in the outer parts. In other studies, the decrease of the hypothetical turbulent component is inferred from the observation that near the core’s center, the linewidth becomes approximately constant at a value marginally larger than the sound speed, a result which is interpreted as the turbulent component becoming smaller than the sound speed (e.g., Goodman et al., 1998; Chen et al., 2018).

However, an alternative interpretation of the decrease of the linewidth towards the central parts of the cores is that the nonthermal motions do not correspond mainly to turbulence, but rather to the infall motions. For prestellar cores, then, the amplitude of the infall motions decreases towards the center, eventually becoming subsonic. Thus, the decrease in the linewidth towards the cores’ centers may correspond to the decreasing infall speed, rather than to a drop in the turbulent component. We plan to further investigate this possibility in a future contribution.

6.6 Development of the hierarchical collapse and observed

fragmentation properties of cores

6.6.1 Independence of time for onset of collapse from the

turbulent Mach number

In Sec. 3.4 we found that the time to reach the $n$th level of fragmentation is independent of the turbulent Mach number. Although surprising, this result reaffirms the notion that the role of the turbulence is only to produce the seeds for collapse, but that the mass of the collapsing objects is determined by the temporal evolution of the mean Jeans mass in the cloud as it contracts gravitationally. This result is consistent with the recent numerical finding by Guszejnov et al. (2018) that the degree of fragmentation of a turbulent cloud depends only on the number of Jeans masses it contains, and not on the Mach number of the turbulence. Finally, it is also consistent with recent observational results that the fragmentation level101010Observationally, the fragmentation level has been defined as the number of fragments –compact sources, of size $\la 5000$ AU and masses of at least $0.5\,M_{\odot}$– within a fixed-size field of view centered on the peak of the clump. This definition is not exactly the same as ours, but refers to the same basic process. observed at scales of 0.1 pc does not appear to correlate well with the turbulent Mach number (e.g., Palau et al., 2015; K. Lee et al., 2015; Beuther et al., 2018), and instead correlates with the density of the parent structure, as discussed in the next subsection.

6.6.2 Dependence of the fragmentation level on the density of

the parent structure.

A natural consequence of the GHC scenario is that the fragmentation level within a given structure should be directly proportional to the Jeans number $N_{\rm J}$ (mass of the structure divided by the Jeans mass) of that structure. In addition, given that $N_{\rm J}\propto\rho^{1/2}$, it is then natural that in the GHC scenario a correlation is expected between the fragmentation level of a structure of a given size and the density of such a structure (at the corresponding size scale). This is fully consistent with a number of observational works reporting a tight relation between the fragmentation level or young stellar objects counts and the density of their parent structures (Gutermuth et al., 2011; Palau et al., 2014, 2015, 2018; K. Lee et al., 2015; Mairs et al., 2016, 2017; Sharma et al., 2016; Hacar et al., 2017, 2018; Alfaro & Román-Zúñiga, 2018; Murillo et al., 2018; Pokhrel et al., 2018; Li et al., 2019; Orkisz et al., 2019; Sokol et al., 2019; Zhang et al., 2019). However, the aforementioned correlation between the fragmentation level (assessed by counting compact fragments) and the density of its parent structure (at the corresponding scale where the fragmentation level is assessed) is only expected if the observed structure is already in an advanced stage after the global collapse started (Fig. 4, right panel). Otherwise, the last fragmentation stage might have not been reached yet, preventing the detection of compact fragments (see Sec. 6.6.3).

6.6.3 Time-delay to reach the $n$th level of fragmentation and

observed lack of fragmentation in massive cold clumps

Also in Sec. 3.4, we found that higher fragmentation levels are reached at later times (cf. Fig. 4, right panel; see also Hoyle 1953). This may explain the apparent lack of fragmentation in some massive, cold (infrared-quiet) clumps. While the masses and scaling of the number of fragments with density has been found to be consistent with thermal Jeans fragmentation in a number of clumps (e.g., Palau et al., 2015; Pokhrel et al., 2018), these clumps are often very active sites of star formation (infrared-bright). However, in other works, focusing on colder infrared-quiet clumps at even earlier evolutionary stages, it is found that the actual masses of the fragments are often significantly larger than the Jeans mass (e.g., Bontemps et al., 2010; Wang et al., 2014; Zhang et al., 2015; Csengeri et al., 2017).

This observation may be interpreted as there existing a “delay” in the appearance of the fragments corresponding to the value of the mean Jeans mass in the clump at the time of measurement. The GHC scenario naturally requires this delay to exist, because the fragments must reach a significant density contrast with respect to their parent structure in order to be detected. Indeed, the typical densities of the fragments observed with interferometers are 2-3 orders of magnitude larger than that of the parent structure observed with a single-dish instrument (e.g., Palau et al., 2013).

Therefore, the fragments are not expected to be detected immediately at the onset of their own collapse (even assuming, as we did in Sec. 3.4, that they start with the “typical” density of the rms fluctuation, $\sim 10\times$ the mean density in the parent structure), but rather after their density has grown by at least one or two orders of magnitude, so that the Jeans mass in the fragment will have decreased by a factor of $\sim 3$–10. As a consequence, the measured Jeans mass will be smaller than that corresponding to the time when the fragments began to contract locally, and the fragments themselves will have only recently begun to sub-fragment themselves.

We conclude that the evolutionary nature of the GHC scenario implies that, by the time when the fragments are sufficiently denser than their parent structure as to be singled-out observationally, their Jeans mass will be smaller than that corresponding to the time when they first began to grow, and that they will only be at an early stage of the next level of fragmentation. Another implication is that the new generation of sub-fragments is not expected to exist at random locations in the parent clump, but rather within the fragments already present. Therefore, these lower-mass sub-fragments are predicted in our scenario to have low density contrasts and to be located within the next-lower-hierarchy fragments. This suggests that searches for the next fragmentation stage should be performed within the fragments already present, at high signal-to-noise and high mass sensitivity.

6.7 Comparison with the “gravoturbulent” scenario

Table 1 presents a summary of how the gravoturbulent and the GHC scenarios deal (or do not) with different properties and features of MCs and their SF activity. For reference, the magnetic support scenario (e.g., Shu et al., 1987; Mouschovias, 1991) is also included, as it actually did provide explanations for some observed features that lack a counterpart in the gravoturbulent scenario.

6.7.1 The difference in the flow regimes

The GHC regime differs strongly from the so-called gravo-turbulent (GT) scenario (e.g., Vázquez-Semadeni et al., 2003a; Mac Low & Klessen, 2004; Ballesteros-Paredes et al., 2007; Hennebelle & Falgarone, 2012; Hopkins, 2012; Klessen & Glover, 2016) in that it is intrinsically evolutionary. In contrast, the GT scenario is stationary. In the latter, the supersonic nonthermal motions observed in MCs correspond to supersonic turbulence driven by some external force (e.g., accretion, supernova explosions, bipolar outflows, etc.). In that scenario, the turbulence plays a dual role in the clouds and their substructures, providing global support against the weight of the structures, while locally producing shocks that generate small-scale density enhancements, in which the local Jeans mass (possibly including turbulent pressure) decreases sufficiently as to become smaller than the fluctuation’s mass, causing its collapse. Therefore, the clouds as a whole are in a nearly stationary state, being in approximate virial equilibrium, and only slowly losing mass to star formation at the rate permitted by the turbulent formation of collapsing small-scale structures (e.g., Vázquez-Semadeni et al., 2003a; Mac Low & Klessen, 2004).

However, numerical simulations have shown that, in general, the turbulence generated by the cloud assembly process (e.g., Koyama & Inutsuka, 2002; Audit & Hennebelle, 2005; Heitsch et al., 2005; Vázquez-Semadeni et al., 2006, 2007, 2010; Wareing et al., 2019) is only moderately supersonic (Mach numbers ${\cal M}_{\rm s}\sim 3$) with respect to the dense, cold gas, in contrast with the highly supersonic regimes (${\cal M}_{\rm s}\ga 10$) indicated by observations of MCs (e.g., Larson, 1981; Heyer & Brunt, 2004). As a consequence, the turbulence self-consistently generated during the cloud assembly process eventually becomes insufficient to support the clouds against their self-gravity (e.g. Vázquez-Semadeni et al., 2007; Heitsch et al., 2009). In addition, Clark & Bonnell (2005) have concluded from analysis of the energy budget of turbulent density flutuations in numerical simulations of decaying turbulence, that the Jeans mass in the fluctuations does not directly become smaller than the fluctuation’s mass. Instead, the fluctuations grow by the self-gravity of the larger-scale region.

In view of the above, and in contrast with GT, in GHC the chaotic, multi-scale infall motions are assumed to dominate at all scales (except, perhaps, in a fraction of the low-column density clouds and clumps that appear strongly super-virial; see Sec. 6.4.1), and the truly turbulent motions (those that are fully disorganized and might in principle provide a turbulent pressure) appear to be fed by the gravitational collapse of the structure in which they are observed, but they do not manage to significantly retard the collapse, even though they may approach levels approaching those of virial balance (Guerrero-Gamboa et al., in prep.); this appears to be a consequnce of the strongly dissipative nature of turbulence, which is usually not included in virial treatments of the turbulent contribution. The only role of these turbulent motions is to provide nonlinear density fluctuations that, when they become unstable due to the global reduction of the Jeans mass caused by the global collapse, begin to collapse on their own. Then the pattern repeats itself inside these new, smaller and denser objects. The mechanism can be considered a joint mass and energy cascade driven by self-gravity.

Such a regime had already been envisioned by Field, Blackman & Keto (2008), although they imagined that the flow was virialized at all scales. Instead, Ballesteros-Paredes et al. (2011a) proposed that the flow is dominated by infall motions at all scales, eliminating the need for the flow to first virialize at each scale and then become unstable again to proceed to the next stage of collapse. The latter authors pointed out that the free-fall velocity is only a factor of $\sqrt{2}$ larger than the virial velocity, and therefore free-falling clouds follow the same scaling relation between velocity dispersion ($\sigma_{v}$), size ($R$) and column density ($\Sigma$) ($\sigma_{v}/R^{1/2}\sim(G\Sigma)^{1/2}$; Heyer et al., 2009) as virialized clouds, within typical observational uncertainties.

Additionally, the fact that any given object only lasts as a coherent unit (before fragmenting again) while it contains one to a few Jeans masses implies that very large velocities are in general not observed within a coherent structure. The large velocities associated with late stages of the collapse should be observable as fragment-to-fragment velocities instead. A similar conclusion has been reached by Hacar et al. (2016) concerning the CO linewidths in nearby clouds through an analysis of the correlation between the linewidths and the line opacity.

6.7.2 Inconsistency between turbulence and generalized equipartition

An important point to note is that, in the GT scenario, turbulence is assumed to be in approximate virial equilibrium with self-gravity at all scales, and the Larson (1981) linewidth-size scaling relation, $\sigma_{v}\propto R^{1/2}$ (or, equivalently, ${\cal L}\approx$ cst.) is assumed to be the result of supersonic turbulence with a Burgers’ energy spectrum of the form $E(k)\propto k^{-2}$ (e.g., McKee & Ostriker, 2007), where $E(k)$ is the specific kinetic energy per wavenumber interval, and $k$ is the wavenumber. However, this assumption disregards the fact that the Larson linewidth-size scaling only holds for objects of similar column density, while, when objects with a wide range of column densities are considered, a Heyer-like relation of the form $\sigma_{v}\propto\sqrt{\Sigma R}$ (or, equivalently, ${\cal L}\propto\Sigma^{1/2}$) holds (Heyer et al., 2009; Ballesteros-Paredes et al., 2011a, 2018; Traficante et al., 2018a). This relation arises from approximate energy equipartition between self-gravity and non-thermal motions, either because of virialization or free-fall (cf. Secs. 3.2 and 6.4). This implies that the column density-dependence of the Larson ratio ${\cal L}$ observed in MCs and their dense substructures is inconsistent with Burgers turbulence when objects of a wide range of column densities are considered, as was further discussed in Sec. 2.2.

6.7.3 The role of accretion: filaments, hub accretion,

development of the massive stellar component, and brown dwarfs

While in GT filaments are assumed to be the result of shocks, with no need for developing a longitudinal flow (e.g., Padoan et al., 2001), in GHC the filaments constitute the very accretion flow from large to small scales, akin to rivers carrying the material down the gravitational potential from high “altitudes” to “lakes” (the stars). Because the filaments may extend over scales of up to tens of parsecs, accretion motions also extend over such large scales, albeit with an extremely anisotropic structure, rather than through the common picture of spherical, homologous collapse.

In addition, if the accretion onto the hubs is Bondi-Hoyle-type with respect to the local environmental conditions (cf. Sec. 6.2.2), then the accretion rate onto the hubs should increase over time as the mass of the hubs themselves increases. This may be at least part of the reason massive stars appear late in the simulations. Higher hub accretion rates are expected to require stronger stellar feedback in order for the filaments to be dispersed, and thus result in the formation of more massive stars within more strongly accreting hubs. This may explain the apparent lack of massive stars in regions where the star formation activity is still rather scattered, even if the total mass content is already large, such as the IRDC M17 SWex (e.g., Povich et al., 2016), the apparent excess of massive stars in regions of extremely high SFRs, such as the Galactic Center cluster and the Arches cluster (e.g., Lu et al., 2013; Hosek et al., 2019), and the observed correlation between the mass of the most massive star in a cluster and the cluster mass (e.g., Weidner & Kroupa, 2006). We shall investigate this mechanism in future contributions.

Finally, it is worth noting that the secondary star formation occurring in the cores along the filaments (rather than the central hubs) must involve lower local accretion rates, and moreover may have their gas supply interrupted by the feedback from nearby main hubs already forming massive stars. Thus, the secondary collapses occurring in the filaments are consistent with suggestions that brown dwarf formation may require secondary formation in the periphery of the main star-forming hubs (e.g., Thies et al., 2015).

6.8 Comparison with other scenarios

The GHC scenario is much closer to the competitive accretion one (e.g., Bonnell et al., 2001; Bonnell & Bate, 2006), which already includes small- (at the single star or small stellar group level) and large-scale (at the level of the accretion flow onto the small star-forming sites) collapse flows. The main difference is that the competitive accretion scenario in general overlooks the evolution of the MCs that occurs over timescales of order several megayears, in particular the global cloud contraction and the evolution of cloud properties and the SFR, since the numerical simulations of that process typically only consider parsec-scale clumps and durations $\la 1$ Myr, because they are often aimed at studying the details of the stellar products (multiplicity, IMF, brown dwarf production, etc.; e.g., Bate, 2012; Kuznetsova et al., 2015), rather than at the evolution of the GMC-scale system, over timescales of a few tens of megayears, from their formation to their dispersal, passing through their gravitational contraction, as done in converging-flow simulations (e.g. Vázquez-Semadeni et al., 2007, 2010; Heitsch & Hartmann, 2008; Hennebelle et al., 2008; Banerjee et al., 2009; Colín et al., 2013; Carroll-Nellenback et al., 2014; Körtgen et al., 2016; Wareing et al., 2017a; Zamora-Avilés et al., 2018). An exception is the study by Smith et al. (2009), which did consider the accretion onto the star-forming clumps and their mass evolution, concluding that the massive stars form when the clumps have become massive enough themselves, similarly to Vázquez-Semadeni et al. (2017). But at the local level, our scenario and simulations are fully consistent with competitive accretion.

Also, as discussed in Gómez & Vázquez-Semadeni (2014) and Sec. 3.1, and shown in Fig. 6, the GHC scenario is consistent with the “conveyor belt” scenario (Longmore et al., 2014), since the filamentary accretion flows serve as funnels for the gas to flow from the large (cloud and clump) scales where the gas is initially distributed to the small (core and hub) scales in which it resides once it has been gravitationally compressed to high densities and forms stars. These filaments accrete from the cloud scale, and onto the cores and hubs. Moreover, secondary star formation is observed to occur in the simulations when they have accreted enough mass that they become locally gravitationally supercritical (see Figs. 2 and 3 of Vázquez-Semadeni et al., 2017).

Finally, the scenario outlined here is fully consistent with the “dynamical model for SF” presented by Elmegreen (2015), in which the ISM is in a dynamical cycle between gravitational contraction followed by feedback-induced dispersal. In that model, the dispersal occurs on the free-fall time of the star-forming gas. In our scenario, this time is longer than the free-fall time due to the anisotropy of the infalling flow, which adopts a filamentary shape, and thus collapses on timescales of the order of the free-fall time times the aspect ratio of the structures (Toalá et al., 2012; Pon et al., 2012). But the overall dynamical state is similar to that outlined by Elmegreen (2015), with the infalling gas having virial parameters $\alpha\sim 1$, while the dispersing gas has $\alpha\gg 1$ (e.g., Fig. 11 of Colín et al., 2013). Observations showing that much of the molecular gas mass has large values of the virial parameter (e.g., Heyer et al., 2001; Kauffmann et al., 2013; Leroy et al., 2015; Heyer & Dame, 2015; Traficante et al., 2018a) are thus consistent with the numerical result that most of the gas in an initially contracting MC is dispersed by stellar feedback before it makes it into stars. This feature thus allows the GHC model to avoid the Zuckerman & Palmer (1974) conundrum that free-falling MCs should have a SFR much larger than observed.

6.9 Deconstructing recent criticisms to GHC

The GHC scenario has recently encountered criticisms, which mostly originate from the original SFR and line-shift-absence conundrums (cf. Sec. 2.1). The modern version of the SFR conundrum is that the star formation rate per free-fall time, $\epsilon_{\rm ff}$, is very low, of the order of $\sim 1\%$, in all molecular structures, from GMCs to massive star-forming clumps (e.g., Krumholz & McKee, 2005; Krumholz & Tan, 2007). The modern version of the line-shift-absence conundrum is the statement that traditional infall signatures, such as blue-excess or inverse P-Cygni line profiles are not observed at large scales in MCs (see, e.g., Sec. 4.6 of Murray & Chang, 2015). In addition, Krumholz et al. (2019, hereafter K19l see also Krumholz & McKee 2019) have recently argued that the kinematic signatures of clusters do not correspond to what would be expected from the GHC scenario. All of these criticisms arise from an incomplete consideration of the full phenomenology of GHC, in particular its evolutionary and extremely non-isotropic nature, as well as the ubiquitous presence of accretion at all scales.

6.9.1 The efficiency per free-fall time

The quantity $\epsilon_{\rm ff}$, which is actually an efficiency rather than a rate, is defined as the fraction of a cloud’s mass that is converted into stars over a free-fall time (Krumholz & McKee, 2005); that is,

[TABLE]

where $\dot{M}_{*}$ is the SFR, $M_{\rm gas}$ is the molecular gas mass, and $\tau_{\rm ff}$ is the free-fall time. In practice, $\dot{M}_{*}$ is sometimes measured as $\dot{M}_{*}\approx M_{*}/\Delta t_{\rm SF}$, where $M_{*}$ is the instantaneous stellar mass in an embedded cluster and $\Delta t_{\rm SF}$ is the duration of the star-formation epoch, often taken as the stellar age spread, typically a few megayears (e.g., Evans et al., 2009; Lada et al., 2010; Povich et al., 2016; Retes-Romero et al., 2017). In this case,

[TABLE]

Other studies estimate the SFR from the IR luminosity (e.g., Vutisalchavakul et al., 2016; Liu et al., 2016a).

In either case, eq. (34) may be suitable over large spatial regions or long timescales that involve a large number of clouds, over which the SFR can be averaged. Indeed, Zamora-Avilés & Vázquez-Semadeni (2014) have performed temporal and ensemble averages of their accelerating-star formation model (discussed here in Sec. 5.4.2), and shown that these averages provide good fits to the observed correlation between SFR and dense gas mass, both at the individual cloud scale (when a single-cloud model is temporally averaged), and at the cloud-ensemble scale (when the model is averaged over time and over a cloud-mass spectrum; see Fig. 5 of that paper).

On the other hand, the definition of $\epsilon_{\rm ff}$ in eq. (34) will give a serious underestimation of the efficiency if the present free-fall time is used to estimate it in an individual, evolving star-forming region, since in this case the collapse rate and the star formation rate are increasing strongly, and the free-fall time is decreasing strongly. This can be clearly seen from eq. (11), from which one can write an evolution equation for the ratio of the instantaneous to the initial free-fall time of a collapsing object, which reads

[TABLE]

where, as in Sec. 3.3, $\tau$ is the time in units of the initial free-fall time ($\tau_{{\rm ff},0}$), and $\tau_{\rm ff}(\tau)$ is the instantaneous value of the free-fall time at time $\tau$. From this expression, the mean free-fall time up to time $\tau$ is given by

[TABLE]

For example, toward the end of the collapse, $\tau\rightarrow 1$ and $\tau_{\rm ff}/\tau_{{\rm ff},0}\rightarrow 0$. However, $\langle\tau_{\rm ff}/\tau_{{\rm ff},0}\rangle_{\tau=1}\approx 0.68$, meaning that the characteristic collapse timescale is not shorter than 2/3 of the initial free-fall time. Using the final free-fall time severely underestimates the mean timescale, spuriously producing very low values of $\epsilon_{\rm ff}$ when it is estimated using eq. (34).

As an example, if a star-forming clump is $\sim 100\times$ denser than the cloud it formed from, its free-fall time will be shorter than the initial one by a factor $\sim 10$, and shorter than the mean one by a factor $\sim 6.8$. Interestingly, this factor is of the order of magnitude of the factor by which K19 suggest that the duration of the star-formation activity is longer than the measured free-fall time (a factor $\sim 10$). The interpretation in the GHC scenario is that the present-day free-fall time is up to $10\times$ or more shorter than the mean one over the development of the cluster up to its present stage. The clouds have been evolving at the free-fall rate, but this rate was much lower over most of the past history of a forming cluster than it is at present. This is illustrated, for example, by the evolutionary track in Fig. 8, which shows that the bursting stage of the model cloud lasts only a small fraction of the whole evolutionary process (less than 1 Myr in that specific example).

It is worth pointing out that this bias is analogous to that produced in the estimate of the expected fragmentation level in massive clumps, as discussed in Sec. 6.6.3. In that case, it was a consequence of the usage of the present-day Jeans mass, rather than the initial one. Here, it is a consequence of using the present-day $\tau_{\rm ff}$ rather than the characteristic one.

It is important to note that the above discussion applies to individual MCs that are observed at random stages of their evolutionary process. At the scale of hundreds of parsecs or more, of the type of regions observed in nearby external galaxy surveys (e.g., Bigiel et al., 2008; Leroy et al., 2008; Onodera et al., 2010; Schruba et al., 2010; Liu et al., 2011), which contain ensembles of GMCs as well as more diffuse gas, the GHC interpretation of the low observed $\epsilon_{\rm ff}$ is different. In this type of studies, the SFR and the gas mass are estimated by means of the region brightness in gas or star-formation tracers, such as CO and Hα emission, respectively. In this case, the reason for the low observed values of $\epsilon_{\rm ff}$ is the same as the solution to the SFR conundrum: since most of the mass in initially contracting MCs is prevented from actually forming stars, and is instead dispersed or evaporated, large regions are expected to contain much of this “molecular debris” from a previous SF episode. Note, however, that this requires that an important fraction of the molecular mass is dispersed rather than evaporated. This can be considered a prediction of the scenario, which should be tested in simulations of the full life cycle of MCs, in the presence of various feedback mechanisms. Existing simulations using photoionising radiation feedback (e.g., Colín et al., 2013; Zamora-Avilés et al., 2018; Haid et al., 2019) suggest that indeed a substantial fraction of dense gas survives.

6.9.2 The absence of infall signatures

The argument is often made that MCs cannot be in global collapse because characteristic infall profiles are generally not observed at the cloud scales ($\ga 10$ pc; e.g., Murray & Chang, 2015, hereafter MC15). The GHC explanation for this is that the collapse flow is far from isotropic, and instead proceeds mostly along filamentary structures. This collapse flow is indeed observed, except not so much as a velocity gradient along the line of sight which produces an infall line profile, but rather as longitudinal flows along the filaments, which feed the central clumps (e.g., Sugitani et al., 2011; Kirk et al., 2013; Fernández-López et al., 2014; Motte et al., 2014; Peretto et al., 2014; Tackenberg et al., 2014; Jiménez-Serra et al., 2014; Hajigholi et al., 2016; Wyrowski et al., 2016; Juárez et al., 2017; Rayner et al., 2017; Lu et al., 2018; Baug et al., 2018; Gong et al., 2018; Ryabukhina et al., 2018; Dutta et al., 2018; Chen et al., 2019). Away from the filaments, the flow is directed mainly towards the filaments themselves rather than towards the central hubs (see the left panel of Fig. 6 and the observational results by Shimajiri et al., 2019). Therefore, in the GHC scenario, the flow at cloud scales is highly chaotic and a simple signature such as typical infall line profiles, which is based on the assumption of a roughly spherical structure and flow, is not generally expected, although, as mentioned in Sec. 2.1, recent studies at the clump (Barnes et al., 2018) and GMC (Schneider et al., 2015) scales have indeed found infall signatures by combining optically thick 12CO and optically thin 13CO lines.

6.9.3 The linewidth-size relation at large and small scales

MC15 have also argued that the different observed exponents $p$ of Larson’s (1981) linewidth-size relation at the scale of whole GMCs ($p\sim 0.5$; e.g., Solomon et al., 1987; Heyer & Brunt, 2004) are larger than at the scale of massive dense cores ($p\sim 0.25$; e.g., Plume et al., 1997; Caselli, & Myers, 1995). MC15 arrive to this conclusion because they build a model for the “adiabatic heating” of turbulence from gravitatioanl contraction, in which they find two different regimes: an external region in which $p\sim 0.5$ and an internal one, dominated by the gravity of a central stellar object, where $p\sim 0.25$ (see also Murray et al., 2017). Thus, MC15 suggest that the lower values of $p$ in the dense cores are indicative of gravitational contraction at the clump scale, while the higher values at the GMC scale are indicative of supersonic turbulence rather than gravitational contraction.

However, Ballesteros-Paredes et al. (2011a) have shown that both GMCs and massive clumps do follow a unique scaling of the form of eq. (1). That is, the value $p\sim 0.5$ applies for roughly-constant-column density objects, such as molecular clouds, but when the column density varies, the linewidth-size relation is shifted, as also shown in numerical simulations by Camacho et al. (2016). Thus, the shallower linewidth-size relation for massive clumps can be understood as a consequence that the larger clumps have lower column densities. That this is so is demonstrated by the fact that, when the clumps and the GMCs are plotted in the ${\cal L}$–$\Sigma$ diagram, they all fall in the same scaling relation.

6.9.4 Radial motions of cluster members

It has been recenty claimed by K19 that the GHC scenario should produce radial motions in the members of stellar clusters formed by the collapsing clouds which, they claim, are not observed. Thus, they conclude that the GHC scenario does not match the observed young cluster kinematics. This line of argument has two problems:

First, K19 start by stating that the GHC scenario should produce radial infall motions in the stellar products during the early stages of the process, and then argue that these inwards radial motions are not observed, thus concluding that collapse cannot be occurring. This claim seems to arise from an implicit assumption of roughly spherical symmetry, which is however not backed up by the very simulations that follow the collapse process (Kuznetsova et al., 2015, 2018a; Vázquez-Semadeni et al., 2017). In the simulations, the structure of the forming clusters is hierarchical (or fractal), having been inherited from their parent clouds. Thus, the clusters consist of groups and sub-groups, and the motions are a combination of the local plus the global collapse flows, on a highly filamentary substrate (see, e.g., Figs. 2 and 3 of Vázquez-Semadeni et al., 2017). Thus, the claim of radial inward motions arises from an oversimplification of the hierarchical collapse mechanism, rather than from the actual kinematics observed in the simulations of GHC.

Second, for more advanced stages, K19 argue that, once the gas is expelled by the stellar feedback, up to 90% of the stars should be moving radially away from the dense cluster, and that thus, the GHC scenario should produce either inward or outward radial motions of the cluster members. Then they quote Ward & Kruijssen (2018), Kuhn et al. (2019), and Kounkel et al. (2018) as examples that Gaia shows no such radial motions in Orion and other complexes, again concluding that the observations contradict the GHC scenario.

This line of argument is actually confusing, for two reasons. The first is that the expansion due to gas expulsion should occur independently of the scenario, since it only involves a shallowing of the gravitational potential well, caused by the feedback, independently of the mechanism that assembled the cluster-forming gas. So, the expansion should be present regardless. The second is that both Kuhn et al. (2019) and Kounkel et al. (2018) do report expansion motions, and in fact the former authors report typical expansion velocities of $\sim 0.5\,{\rm~{}km~{}s}^{-1}$, and even a radial expansion velocity gradient, so the reference to these papers by K19 against expansion seems contradictory. In any case, the simulations with feedback presented by Vázquez-Semadeni et al. (2017) also show expansion (again, see Figs. 2 and 3 of that paper), although a detailed comparison with the observations remains to be done.

Concerning Ward & Kruijssen (2018), these authors do claim that no expansion is seen in their data. However, they consider much larger regions, several tens of parsecs across. As mentioned in Sec. 5.5, at size scales this large, the coherence of the infall motions is lost because the infall timescale is larger than the cloud-destruction timescale from the feedback, so no clear signature of the infall is to be expected. This is aided by the multi-scale and filamentary nature of the infall motions, which contribute to the non-appearance of a single-focused collapse at those scales. On the other hand, Ward & Kruijssen (2018) do conclude that their observations support a hierarchical star formation model.

We thus conclude that the claim by K19, that observed cluster kinematics do not support the GHC scenario, is unfounded.

6.10 Caveats

In this paper we have presented several calculations of the evolution of clouds as they undergo gravitational contraction, but of course, these calculations are only approximate, and are presented mostly as proof of concept. The main approximations we have used, and which must be improved upon in more detailed calculations (such as those performed semi-analytically in the ZV14 model) are as follows:

First, we have assumed spherical symmetry and uniform density by using eq. (10) for the evolution of the radius of a collapsing object. In reality, clouds are expected to be flattened, and thus to evolve more slowly than a spherical cloud with the same density (Burkert & Hartmann, 2004; Toalá et al., 2012; Pon et al., 2012; Zamora-Avilés, Vázquez-Semadeni & Colín, 2012). Moreover, due to the presence of fluctuations, the clouds contract non-homologously, as shown by both similarity studies and numerical simulations (e.g., Larson, 1969; Penston, 1969; Hunter, 1977; Shu, 1977; Whitworth & Summers, 1985; Foster & Chevalier, 1993; Mohammadpour & Stahler, 2013; Keto et al., 2015; Naranjo-Romero et al., 2015). Thus, we have substituted the calculation of a single, non-homologous contraction, by the consideration of a sequence of uniform spheres of different densities that collapse on different timescales.

Second, because of this treatment, we have neglected accretion onto an object at a given scale from its parent object at a larger scale. Accretion from the next-larger scale has been shown by numerical simulations to be an essential part of the multi-scale collapse process (e.g., see Fig. 7 in Vázquez-Semadeni et al., 2009, and Sec. 4.3.2 of Vázquez-Semadeni et al. 2010). Also, the clumps lose mass by conversion to stars and by the effect of stellar feedback. These important aspects of the evolution remain to be included in a more thorough analytical treatment. Numerical simulations show that clumps defined by density thresholds rather than by fixed mass (as assumed in Sec. 3.2) evolve along different trajectories than those shown in Fig. 2, turning back to lower column densities as they begin to lose mass (Camacho et al. 2019, in prep.).

Third, we have neglected deviations from spherical or circular geometries as the collapse proceeds. As proposed by Gómez & Vázquez-Semadeni (2014) and discussed in Sec. 5.3, the infall of the large scales develops filamentary flows, which feed the faster-evolving, smaller-scale, more roundish objects. However, this should not introduce significant errors in the calculation of the collapse timescales since, as discussed in Toalá et al. (2012), the collapse timescale of a flattened or filamentary structure is the same as that of a spherical structure containing the same mass over size scale equal to the largest spatial extent of the sheet or filament. Since the filaments basically form from the collapse of the large-scale cloud, their timescale will be the same of the initial, closer-to-spherical diffuse parent cloud.

Finally, we have neglected any effects of support provided by the Galactic rotation and shear (i.e., the Toomre criterion). This implies that our discussion is limited to scales where the shear induced by the Galactic rotation is negligible, i.e., a few hundred parsecs.

7 Summary and conclusions

7.1 Summary

In this paper we have presented a complete and unified description of the GHC model, starting from the formation of MCs as thin, cold atomic sheets, through their growth, evolution of their SFR, until their destruction by the massive stars they form, preventing most of their mass from being converted to stars, and thus keeping the overall SFE low. The model describes the phenomenology observed in numerical simulations of MC formation in the warm atomic medium. Figure 13 shows a schematic representation of the evolution of a $10^{4}$–$10^{5}M_{\odot}$ GMC for typical Solar neighborhood conditions, assuming it starts from WNM gas exclusively. The next subsections briefly summarize the evolutionary sequence.

7.1.1 Cloud formation and the onset of collapse

In this scenario, MCs at galactocentric distances where the majority of the volume of the galactic disk is occupied by diffuse gas ($R\ga 6.5$ kpc; Koda et al., 2016) initiate their life cycle by compressions in the predominantly diffuse medium, which nonlinearly trigger a transition to the cold atomic phase (Hennebelle & Pérault, 1999; Koyama & Inutsuka, 2002). Clouds formed by this type of compression tend to be planar, because the collision of streams or shock fronts are most often two-dimensional surfaces, which grow by accretion of gas from the warm phase. The clouds are bounded by a transition front, where the accreted gas undergoes the transition from the warm to the cold phase (Vázquez-Semadeni et al., 2006; Banerjee et al., 2009).

Because the cold phase is $\sim 100\times$ denser and colder than the initial warm phase, the Jeans mass in this material is $\sim 10^{4}\times$ smaller than in the warm phase, and thus the cold gas quickly (in a few to several megayears, depending on the coherence scale of the converging motions) acquires masses much larger than its Jeans mass (Gómez & Vázquez-Semadeni, 2014). The accumulation of gas produces turbulence, but generally it is only moderately supersonic (Vishniac, 1994; Walder & Folini, 2000; Koyama & Inutsuka, 2002; Heitsch et al., 2005; Vázquez-Semadeni et al., 2006; Wareing et al., 2019, ${\cal M}_{\rm s}\sim$ a few), implying that in general it is insufficient to prevent the cloud from quickly engaging in global gravitational collapse, at which point the nonthermal kinetic energy becomes “locked” to the gravitational energy, but is never capable of retarding the collapse (Vázquez-Semadeni et al., 2007).

7.1.2 Fragmentation and filament formation

Because the cloud acquires masses much larger than the Jeans mass, the collapse becomes nearly pressureless, and thus it proceeds fastest along the shortest dimension of the cloud, reducing its dimensionality (Lin et al., 1965), leading to the formation of filaments (Heitsch et al., 2008b; Gómez & Vázquez-Semadeni, 2014). Within the filaments, roundish density fluctuations grow faster than the rest of the filament because the filaments have longer collapse timescales (Toalá et al., 2012; Pon et al., 2012). Thus, these clumps produce a deeper gravitational potential, and the rest of the filament begins to accrete onto the clump, developing a longitudinal flow along the filament and onto the clumps (Gómez & Vázquez-Semadeni, 2014). It is important to remark that no strong shock is produced at the axis of the filament because a shock only develops in a gravitationally contracting structure when the density develops a singularity, but this is avoided by the longitudinal mass flux along the filament towards the clumps, which apparently leads to the establishment of a near-stationary state in the filaments (Naranjo-Romero et al.  in prep.). The absence of a central shock is indeed observed in numerical simulations of gravitationally-formed filaments (Fig. 6), and may constitute a diagnostic to distinguish between gravity- and turbulence-formed filaments.

As the mean density in the cloud increases, the mean Jeans mass decreases (eq. [14]). Assuming that the density PDF of the turbulent fluctuations (not of the collapsing objects) retains its lognormal shape and its mean Mach number, the density PDF evolves simply by shifting to higher densities while retaining its shape (Zamora-Avilés, Vázquez-Semadeni & Colín, 2012). Moreover, because the clouds are sheet-like and the turbulent density fluctuations are nonlinear, Hoyle-like fragmentation can occur without concern for the objections raised by Tohline (1980), as discussed in Sec. 3.1.

7.1.3 Increase of the SFR and the mean stellar mass formed

At the typical density of the turbulent density fluctuations, less massive fluctuations become unstable at later times (Fig. 3). However, the first fragments that begin to collapse correspond to the most extreme density fluctuations that have a mass larger than their corresponding local Jeans mass. Therefore, the first fragments to collapse have the lowest masses. Subsequently, as the mean density increases, the local collapses can occur for densities less distant from the mean. Since there is more mass at these densities, larger-mass fluctuations can collapse, and so the typical mass of the collapsing objects increases with time. Simultaneously, the SFR increases with time (Zamora-Avilés, Vázquez-Semadeni & Colín, 2012; Zamora-Avilés & Vázquez-Semadeni, 2014; Vázquez-Semadeni et al., 2018), and so the star-forming clouds evolve towards larger SFRs while increasing the mean mass of the stars they form (see Fig. 10 here, taken from Vázquez-Semadeni et al. (2017)). This implies that there should be a correlation between the most massive star in a star-forming core or clump, and the clump’s mass. In turn, this implies the existence of a correlation between the mass of the most massive star in a cluster and the cluster mass, which is indeed observed (e.g., Weidner & Kroupa, 2006).

7.1.4 Energy balance evolution in collapsing clumps

Moreover, as the clump evolves, so does its virial parameter. According to the discussion in Ballesteros-Paredes et al. (2018) and Sec. 3.2, the virial parameter contains an inertial, possibly compressive, component111111Recall that this component need not consist of turbulence that supports the cloud. It can just as well consist of the external converging flows that produce the clump. and a gravitationally-driven infall one ($v_{\rm g}$). Because the collapse of fragments starts at an initial radius $R_{0}$, the infall speed is given by eq. (5), implying that the gravitational component is expected to be subvirial over a significant fraction of the clump’s evolution. Thus, depending on the initial value of the ratio of the inertial to the gravitational components, the core may appear super- or sub-virial.

Moreover, at a given size scale, the more massive (i.e., denser) cores will have a larger value of $v_{\rm g}$, and thus their total linewidth will tend to be dominated by this component, which may nevertheless be sub-virial. This may explain the often-observed trend for more massive cores to be more sub-virial (e.g., Kauffmann et al., 2013; Ohashi et al., 2016; Sanhueza et al., 2017; Traficante et al., 2018b). We will explore this scenario in a future publication.

7.1.5 Local cloud disruption

As massive stars begin to form in the region, they begin to feed back sufficiently strongly on their parent clump. This has the effect of first disrupting the filaments feeding the central star-forming core or hub, and later the hub itself (Fig. 11) or, alternatively, allowing the gas in the hub to be exhausted because the external gas supply has been interrupted. This process leads to a gradual reduction, and possibly eventual full local supression of the SFR. Numerical simulations of the process show that the final SFEs are $\sim 30\%$ (e.g., Colín et al., 2013; Iffrig & Hennebelle, 2017). Such values are reasonable, considering that they refer to the final value of the SFE, for which there is no observational determination, because by the time a cluster appears devoid of gas, it is not possible to observationally know the total amount of gas that went into its formation. Currently, the issue that still needs to be properly assesed is the fact that the effect of the feedback depends strongly on the type of feedback and the precise location of the injection with respect to the clouds (e.g., Dale, 2015; Iffrig & Hennebelle, 2015).

7.2 Resolution of the old conundrums

A crucial question that may be asked of the GHC scenario is whether it provides an answer to the decades-old conundrums that originally dispelled the scenario that star-forming molecular clouds may be in global gravitational collapse (Liszt et al., 1974; Goldreich & Kwan, 1974). But indeed it does. The SFR conundrum (Zuckerman & Palmer, 1974, cf. Sec. 2.1 here), that the SFR should be much larger than the observed value if the clouds are in global gravitational collapse, is avoided by the fact that the star-forming regions are destroyed by the first very massive stars that form (larger numbers of massive stars may be necessary for destroying more massive clouds; Franco, Shore & Tenorio-Tagle, 1994), so that the majority of the cloud mass never manages to be incorporated into stars, but is rather returned to the ambient medium. Indeed, numerical simulations of clouds up to $\sim 10^{5}M_{\odot}$ including stellar feedback (ionizing radiation, winds and SNe) show that the clouds are readily destroyed, with final SFEs that depend on the type and location of feedback considered (e.g., Dale et al., 2012, 2013; Krumholz & Thompson, 2012, 2013; Colín et al., 2013; Iffrig & Hennebelle, 2015, 2017; Körtgen et al., 2016; Wareing et al., 2017a, b; Seifried et al., 2018). Note that the typical quoted efficiencies for GMCs, of a few percent, refer to the SFE of objects that are still dominated by the gaseous mass, so by definition they must be small.

The other conundrum, which we referred to as the line-shift conundrum (Zuckerman & Evans, 1974), is avoided in the GHC scenario because such line offsets are not expected to be observed in denser gas, since this gas acquires a filamentary morphology, and so there is no roughly spherical envelope around the collapse centers (the hubs). Instead, the gas flows from the diffuse regions goes onto the filaments, and then longitudinally along the filaments onto the hubs (see Fig. 6). In addition, such shifts are beginning to be marginally observed in more diffuse molecular gas (12CO and 13CO) (Barnes et al., 2018).

It should be emphasized that many of the objections against the GHC scenario arise from a failure to appreciate the self-similar and non-homologous nature of a realistic gravitational collapse, which implies that regions farther away from the center can be considered to be at an earlier evolutionary stage than the more central ones. This is best illustrated by a similarity description of spherical gravitational collapse (e.g., Larson, 1969; Penston, 1969; Shu, 1977; Hunter, 1977; Whitworth & Summers, 1985), in which the spatial and temporal variables are merged into a single similarity variable, defined by $\xi\equiv r/c_{\rm s}t$; that is, the radial position is normalized by a time-dependent factor amounting to the distance traveled by a sound wave up to time $t$. Thus, large distances from the collapse center at a given time behave like shorter distances at an earlier time, suggesting that a non-homologously-collapsing spherical core can be thought of as a collection of concentric collapsing shells, each at a different evolutionary stage.

In the analytic solutions for the spherical problem, the infall speed at large radii during the prestellar stage approaches a constant (e.g., Larson, 1969; Whitworth & Summers, 1985). However, in numerical simulations, because the collapse starts from a local fluctuation, the size of the infalling region is observed to grow over time, similarly to the expanding rarefaction front of Shu’s (1977) inside-out collapse (see, e.g., the animations corresponding to Fig. 3 of Gómez & Vázquez-Semadeni, 2014, noting how the longitudinal motions along the filaments progressively extend to ever larger distances). Note, however, that this front lies much further out than Shu’s front, because the infall motions begin much earlier (at the onset of collapse), rather than at the time of the formation of a singularity, as he assumed. So, the internal regions of a finite collapsing prestellar structure exhibit an oustide-in nature (cf. Sec. 4.2), while the external ones exhibit an inside-out one.

Finally, when the feedback from the massive stellar products begins to disrupt the cloud, regions sufficiently distant that the infall-motion front has not reached them yet will be effectively decoupled from the collapse by the feedback, and this material does not participate in the star formation episode of the region, maintaining the global SFR low. Thus, even if clouds are dominated by gravity at large scales, the motions do not necessarily reflect a clear collapse.

7.3 Implications of multi-scale collapse for the gas flow

The condition that MCs and their substructures are undergoing global hierarchical gravitational contraction has the important implication that the gas is flowing from the low- to the high-density regions. This means that the gas changes location as it proceeds to higher densities, and this process takes time. Simply stated, before a hydrogen atom reaches the interior of a star, it must traverse a path where the environmental density increases from the mean density of the Universe to that in the star’s interior.

This apparently trivial consideration has strong implications on notions such as the one that low-mass dense cores may be Jeans-stable and hydrostatic, requiring an external pressure to remain confined. Instead, numerical simulations of spherical collapse on top of a globally Jeans-unstable nearly uniform background (Naranjo-Romero et al., 2015) show that the background is also infalling (accreting) onto the core, and thus the observational truncation process may artificially discard part of the infalling material that provides a ram, rather than thermal, pressure, so that the whole system is undergoing gravitational contraction, rather than in equilibrium.

7.4 Final considerations

The GHC scenario provides a unified framework that allows understanding a wide range of MC properties, such as the apparent virialization at all scales, the formation of filaments, the acceleration of star formation, the low overall SFE, and the Bonnor-Ebert-like structure of cores while allowing them to be contracting structures. It is based on the premise that, opposite to the standard beliefs, star-forming MCs, rather than being near-equilibrium structures, are undergoing nearly pressureless gravitational contraction, albeit in an extremely non-homologous and hierarchical fashion. This is made possible by means of a Hoyle-like fragmentation in the nearly isothermal gas that makes up the clouds, which in turn can occur thanks to the moderate turbulence in the clouds that allows the formation of nonlinear density fluctuations that have shorter free-fall times than the whole cloud. Although they do not begin their collapse until the global collapse has sufficiently reduced the mean Jeans mass in the cloud, they do terminate their own collapse earlier than the whole cloud, allowing for its destruction before most of its mass has made it into stars.

In this paper we have shown the ability of the GHC scenario to correctly describe, at both the qualitative (through the analytical approximation discussed here) and quantitative (through the numerical simulations) the spatial structure of observed clouds, star-forming regions and young clusters, and to make predictions about the evolutionary processes they undergo. It is noteworthy that the GHC scenario bears a strong resemblance to the formation of structure in the early Universe, by means of a cosmic web of filaments in which galaxies form (e.g., Bond et al., 1996; Cantalupo et al., 2014), and within which, in turn, the first clusters form again through filamentary accretion (e.g., Safranek-Shrader et al., 2016).

Claims that the scenario does not agree with observations (e.g., Krumholz et al., 2019; Krumholz, & McKee, 2019) appear unfounded, and to originate from a perceived oversimplification of the mechanisms at play. Observational works performing detailed comparisons at the structural level of clouds and clusters have also reported qualitative and quantitative agreement (e.g., Motte et al., 2018; Getman et al., 2018, 2019; Chen et al., 2019) with the GHC scenario. In future works we plan to continue testing the scenario in a variety of situations.

acknowledgements

We thankfully acknowldge useful and constructive comments from Peter Barnes, Andi Burkert, Bruce Elmegreen, Neal Evans, Roberto Galván-Madrid, Lee Hartmann, Pavel Kroupa, Mark Krumholz, Milos̆ Milosavljević, Nicola Schneider, Alessio Traficante, Chris Wareing, William Wall, and Long Wang. This work has received financial support from CONACYT grants 102488 and 255295 to EV-S, from UNAM-PAPIIT grant IN113119 to A.P., and UNAM-PAPIIT grant IN100916 to G.G.