Supernova Driving. IV. The Star Formation Rate of Molecular Clouds
SSubmitted to ApJ, April 4, 2017
Preprint typeset using L A TEX style emulateapj v. 12/16/11
SUPERNOVA DRIVING. IV. THE STAR FORMATION RATE OF MOLECULAR CLOUDS
Paolo Padoan,
Institut de Ci`encies del Cosmos, Universitat de Barcelona, IEEC-UB, Mart´ı i Franqu`es 1, E08028 Barcelona, Spain; [email protected], Pg. Llu´ıs Companys 23, 08010 Barcelona, Spain
Troels Haugbølle,
Centre for Star and Planet Formation, Niels Bohr Institute and Natural History Museum of Denmark, University of Copenhagen, ØsterVoldgade 5-7, DK-1350 Copenhagen K, Denmark; [email protected] ˚Ake Nordlund,
Centre for Star and Planet Formation, Niels Bohr Institute and Natural History Museum of Denmark, University of Copenhagen, ØsterVoldgade 5-7, DK-1350 Copenhagen K, Denmark; [email protected]
Søren Frimann
Institut de Ci`encies del Cosmos, Universitat de Barcelona, IEEC-UB, Mart´ı i Franqu`es 1, E08028 Barcelona, Spain; [email protected]
Submitted to ApJ, April 4, 2017
ABSTRACTWe compute the star formation rate (SFR) in molecular clouds (MCs) that originate ab initio in anew, higher-resolution simulation of supernova-driven turbulence. Because of the large number of well-resolved clouds with self-consistent boundary and initial conditions, we obtain a large range of cloudphysical parameters with realistic statistical distributions, an unprecedented sample of star-formingregions to test SFR models and to interpret observational surveys. We confirm the dependence ofthe SFR per free-fall time,
SF R ff , on the virial parameter, α vir , found in previous simulations, andcompare a revised version of our turbulent fragmentation model with the numerical results. Thedependences on Mach number, M , gas to magnetic pressure ratio, β , and compressive to solenoidalpower ratio, χ at fixed α vir are not well constrained, because of random scatter due to time andcloud-to-cloud variations in SF R ff . We find that SF R ff in MCs can take any value in the range0 ≤ SF R ff (cid:46) .
2, and its probability distribution peaks at a value
SF R ff ≈ . SF R ff and the scatter in the SF R ff – α vir relation are consistent withrecent measurements in nearby MCs and in clouds near the Galactic center. Although not explicitlymodeled by the theory, the scatter is consistent with the physical assumptions of our revised modeland may also result in part from a lack of statistical equilibrium of the turbulence, due to the transientnature of MCs. Subject headings:
ISM: kinematics and dynamics – MHD – stars: formation – turbulence INTRODUCTION
Star formation in galaxies is ultimately regulated bythe external supply of gas that can rapidly cool and set-tle at the high densities of star-forming regions. The starformation rate (SFR) probed by the Kennicut-Schmidtrelation (Kennicutt 1998), or by the Madau plot (Madauet al. 1998), is essentially determined by the cosmologi-cal environment controlling the mass infall from the cos-mic web onto the galaxies (e.g. Dekel et al. 2009; Forbeset al. 2014) and the radial transport in the galactic diskfrom low density gas at large radii to high density gas atsmaller radii. The dark matter component of this bound-ary condition of galaxies is well understood in the stan-dard ΛCDM model, while the structure of the baryoniccomponent of the infall has been studied with large dy-namic range simulations including gas dynamics. How-ever, galaxies can achieve a global SFR in approximatebalance with the cosmological infall through a varietyof paths, leading to or requiring different galaxy mor-phologies and different dynamical and chemical evolu-tions. For example, the star formation may be smoothly distributed in space or time, or may occur in rapid burstswithin massive clouds. The detailed mode of star forma-tion is only partly dependent on the large-scale boundaryconditions; it is primarily controlled by disk dynamics(e.g. gravitational instability and radial transport) andby the physics of star formation, meaning the specificprocesses responsible for the evolution and fragmentationof cold interstellar medium (ISM) clouds. The theoreti-cal modeling of such processes should result in a physicallaw that predicts the SFR as a function of star-formingcloud parameters. The formation and evolution of galax-ies cannot be fully understood until such a universal SFRlaw is revealed.Despite recent progress in modeling the effect of su-personic turbulence on the fragmentation of star-formingclouds, the development of large dynamic range simula-tions of star formation, and the ever growing number andsize of galactic and extragalactic surveys of star formingregions, SFR laws remain very hard to test. The SFRin MCs is difficult to measure, with different methodsoften yielding very different values. The cloud physical a r X i v : . [ a s t r o - ph . GA ] M a r Padoan et al.parameters on which the theoretical models depend (forexample the ratio of gas and magnetic pressures, or theratio of compressive to solenoidal power of the velocityfield) are difficult to measure as well. Furthermore, re-cent compilations of MC properties including their SFRhave shown that the SFR has only a weak dependence oncloud parameters and a very large scatter (e.g. Murray2011; Evans et al. 2014; Vutisalchavakul et al. 2016), inapparent contradiction of most theoretical models. Be-sides the very large scatter, the observed SFRs also tendto be lower than both numerical estimates and theoreti-cal predictions.On the other hand, idealized numerical simulations ofsupersonic turbulence have been used to test theoreticalSFR laws with some success (Padoan & Nordlund 2011;Federrath & Klessen 2012; Padoan et al. 2012). Theyhave also yielded best-fit values of the free parameters ofthe theoretical models. These simulations tend to givelarger SFR than the observations, even when the turbu-lence is well resolved and relatively strong magnetic fieldsare included. More importantly, they do not yield thelarge scatter in SFR, for fixed physical parameters, foundin the observations. Thus, while our theoretical under-standing of the star formation process has improved sig-nificantly with the interplay of theory and simulations,our ability to reproduce the observations remains limited,casting doubts upon the theoretical scenario of turbulentfragmentation. Nevertheless, we believe the discrepan-cies with the observations are mainly the consequence ofthe limitations of the simulations. More realistic sim-ulations should reproduce the observations and guide arevision of the theoretical SFR models, while retainingthe key idea that star-forming clouds are fragmented bysupersonic turbulence.The main limitation of the idealized simulations usedto test the SFR models is that they capture the funda-mental physics of supersonic MHD turbulence, withoutdescribing a full cloud, its realistic boundary conditionsand its realistic evolution including the cloud formationand dispersion processes. The simulations adopt periodicboundary conditions, so they are interpreted as a charac-teristic piece of a MC, and they are started from idealizedinitial conditions, so they must be evolved for several dy-namical times to pursue a statistical steady state, beforegravity is introduced. Because of the need to first developthe turbulence, gravity has to be suddenly included at alater time, so even the initial evolution under the effect ofgravity is not entirely realistic; the simulation has to berun for at least one free-fall time with gravity, before theSFR is measured. The SFR is initially very low, whengravity is first included, and then gradually increases, inthe best case (but not always) becoming approximatelyconstant for at least a few free-fall times. It is this latetime, approximately (or ‘hopefully’) constant SFR that isusually measured in the simulations, to avoid the initialtransient phase that reflects the imprint of the numeri-cal setup. Thus, the SFR derived in this way from thesimulations does not reflect the initial time variations,because it does not record the low SFR values of theinitial transient phase.Real MCs also have an initial phase when they arefirst assembled and their SFR may be very small. Theymay also become more quiescent after exhausting mostof their dense gas, or reduce their SFR rapidly during a period of cloud dispersion, towards the end of their life-time. But this time evolution of MCs and variability oftheir SFR cannot be described with the idealized simu-lations mentioned above. Thus, it is to be expected thatonly the highest SFR values of MCs are reproduced bythe simulations, while lower values are neglected becausetheir frequency cannot be estimated with the idealizednumerical setups. Furthermore, the physical parametersthat controls the SFR, such as Mach number, virial pa-rameter and gas to magnetic pressure ratio, are usuallyaveraged over the whole computational domain, which isnot the same as measuring them for a whole MC.The alternative approach of simulating a whole cloudfrom the moment of its formation must rely on ad hocand idealized initial conditions, where velocity, densityand magnetic fields can only mimic and grossly misrep-resent the turbulence in real MCs. Because gravity ispresent from the beginning, the star formation processstarts rapidly, and the artificial initial conditions leave astrong imprint on star formation. SFR values derived inthis way cannot be used to faithfully sample the varia-tions of the SFR in real MCs.The only way to model the range of variations of theSFR in MCs is to adopt a numerical setup where MCsare formed ab initio , meaning without ad hoc initial orboundary conditions. This can only be achieved by simu-lating a volume much larger than the size of a single MC(e.g. Dobbs & Pringle 2013; Bonnell et al. 2013; Dobbs2015), which is expensive, but well worth the cost, be-cause a large volume yields a large population of MCs,all emerging from initial and boundary conditions withrealistic probability distributions. With a single largescale run, the SFR evolution can be followed over timein each cloud, generating realistic variations in time andfrom cloud to cloud as well. Very low SFR can then berealized as well, and both average values and scatter ofthe SFR as a function of physical parameters can be com-pared with the observations and with theoretical models,which is the goal of the present work.Such a large scale simulation was already presented inour previous three papers of this series (Padoan et al.2016b; Pan et al. 2016; Padoan et al. 2016a, Paper I,Paper II and Paper III hereafter), where we studied theproperties of the ISM turbulence driven by supernova(SN) explosions. The properties of MCs formed self-consistently in the turbulent ISM of the simulation werefound to agree with those of real MCs from the COFCRAO Outer Galaxy Survey (Heyer et al. 1998, 2001).In this work, we adopt the same numerical setup as a‘laboratory’ to study star formation in MCs, by increas-ing the spatial resolution of the simulation and includ-ing accreting sink particles to describe the formation ofmassive stars. In this new simulation, the turbulenceis driven by the explosion of massive stars whose for-mation is resolved, so both the timing and position ofthe SNe, hence their effect on MCs (Iffrig & Hennebelle2015), is now much more realistic than in the previoussimulation (or any simulation of SN-driven turbulenceto date), where the SNe where generated randomly inspace and time. The main limitation of this setup is thelack of even larger scales, as our computational domainis a cubic volume of 250 pc size with periodic boundaryconditions (hence fixed total mass), neglecting differen-tial rotation and the vertical gravitational potential ofupernova Driving. IV. The Star Formation Rate of MCs 3
Fig. 1.—
Square root of the projected density of the whole 250 pccomputational domain, integrated along the direction of the meanmagnetic field and computed at t = 70 . t = 74 . (cid:12) pc − andsaturates at a maximum column density of 149 M (cid:12) pc − . The pro-jected density is computed using only cells with a number densityabove the threshold n H , min = 200 cm − , to illustrate the massdistribution corresponding to our lower-density MC catalog. the galactic disk. However, such limitation is unlikely toaffect directly the SFR within individual MCs.As in previous works, we normalize the SFR by theratio of the total mass of the system, M , and the free-fall time of its mean density ρ , t ff = (3 π/ Gρ ) / , toobtain a non-dimensional SFR referred to as “SFR perfree-fall time”, SF R ff ≡ − dM/dt/ ( M/t ff ), first intro- Fig. 2.—
The same projected-density fields as in Figure 1, butwith a compressed grayscale intensity range, making the saturationvalue of 149 M (cid:12) pc − darker, thus allowing the sink particles tostand out as dots of the highest intensity value. All sink particleswith mass larger than 1 M (cid:12) are shown (the largest mass is 117.1M (cid:12) ), which amount to 1408 sinks in the upper panel and 2888 sinksin the lower panel. One can clearly see (young) sinks inside thedensest filaments, as well as many (older) sinks that have alreadyleft their parent clouds, including clusters that have cleared theirsurrounding gas thanks to SN explosions of their most massivemembers. duced by Krumholz & McKee (2005). In Padoan et al.(2012), using a large set of AMR simulations, we de-rived a SFR per free-fall time that depended only on thevirial parameter, α vir (the ratio of kinetic turbulent en-ergy and gravitational energy –see definition in Section3), SF R ff , P12 = 0 . − . α / ) (see Section 5.3).Because both the virial parameter and the SFR were de- Padoan et al.rived as global values for the whole computational box,this result is suitable for the development of subgrid mod-els of stellar feedbacks in galaxy formation simulations(e.g. Semenov et al. 2016). In this work, we focus on thedependence of the SFR on the properties on individualclouds, in order to test the theoretical predictions overa large statistical sample and to compare with observa-tional estimates of the SFR in real MCs. By applyinga revision of our turbulent fragmentation model to thephysical parameters of the clouds extracted from the sim-ulation, we derive a new empirical relation between theSFR per free-fall time and the virial parameter that ap-plies to individual clouds, SF R ff ,α = 0 . − . α / , e )(see Section 5.3), where α vir , e is the effective virial pa-rameter of individual clouds (see Section 3), rather thanits estimated value over the whole computational do-main. For reference, we show this new empirical lawas a dashed line in several figures throughout the paper(figures 5, 10, 11, 15), even preceding its derivation inSection 5.3.The paper is organized as follows. In the Section 2, wedescribe the numerical setup, and in Section 3 we defineand compute the physical parameters of MCs selectedfrom the simulation. The cloud SFR is studied in Sec-tion 4, where both average values and scatter are foundto be realistic. We then summarize the turbulent frag-mentation model of the SFR and propose a slight revi-sion in Section 5. In the same section, we also apply ourrevised SFR model to the clouds from the simulation.On average, we find good agreement between the SFRfrom the revised model and the SFR measured directlyfrom the simulation, while the scatter in the model is toosmall. We discuss the origin of the scatter arguing thatit is actually predicted, though not accounted for, by themodel. In Section 6 we compare our results with SFRestimates in real MCs, and in Section 7 we summarizeour work and list the most important conclusions. SIMULATION
This work is based on a new MHD simulation of SN-driven turbulence in a large ISM volume, carried out withthe Ramses AMR code (Teyssier 2002; Fromang et al.2006; Teyssier 2007). The numerical method and setupare discussed extensively in Paper I and are only brieflysummarized here. We simulate a cubic region of size L box = 250 pc, with periodic boundary conditions and atotal mass of M box = 1 . × M (cid:12) . The initial conditionfor this simulation is the final snapshot of our previousSN-driven simulation (see Paper I) before the introduc-tion of self-gravity, at t = 45 Myr. That simulation wasstarted with zero velocity, a uniform density n H , = 5cm − , a uniform temperature T = 10 K, and a uni-form magnetic field B = 4 . µ G, later amplified by theturbulence to an rms value of 7.2 µ G and an average of | B | of 6.0 µ G, consistent with the value of 6 . ± . µ Gderived from the ‘Millennium Arecibo 21-cm Absorption-Line Survey’ by Heiles & Troland (2005). SN explosionswere randomly distributed in space and time (see discus-sion in Paper I in support of this choice), with a rate of6.25 SNe Myr − . The resolution was dx = 0 .
24 pc, witha 128 root grid and three AMR levels. Details about thenumerical setup and the implementation of random SNdriving, tracer particles, and parametrized heating and Fig. 3.—
Star formation efficiency versus time for the wholecomputational volume (black thick solid line), and for individualMCs from the simulation, selected with n H , min = 200 cm − (bluethin lines) and n H , min = 400 cm − (red thin lines). The thinline of each cloud covers 1.68 Myr, the time interval during whichthe clouds are followed (notice that the SFE for five MCs with M cl > × M (cid:12) ( n H , min = 200 cm − ) and for six MCs with M cl > M (cid:12) ( n H , min = 400 cm − ) has been plotted with thickerlines). The vertical dotted lines mark the times of the 10 simulationsnapshots where the MCs were selected, at intervals of 1.5 Myr.The short-dashed and long-dashed lines show SF E versus timefor depletion times t dep = 1 . t = 61 Myr and t dep = 0 . t = 67 Myr, respectively. The twogreen shaded areas show two examples of SF E versus time for t dep = 0 .
05 Gyr, characteristic of the MCs from the simulation,starting at the first and fourth cloud selection times. This clouddepletion time is evaluated as 1 / (cid:104) ( SF R ff /t ff ) (cid:105) averaged over allclouds with non-zero SFR. cooling can be found in Paper I. In the following, we fo-cus on the description of the new features of the currentsimulation.Because we wish to resolve the formation of individualmassive stars, we increase significantly the resolution rel-ative to our previous run. As we restart at t = 45 Myr,we continue to run without self-gravity, and increase theroot grid size to 512 cells, besides adding four AMR lev-els to reach a minimum cell size of dx = 0 .
03 pc. We alsoinitialize 250 million passively advected tracer particles,each representing a fluid element with a characteristicmass of approximately 0.008 M (cid:12) . The tracers record allthe hydrodynamic variables, and are tagged once they ac-crete onto a sink particle, so the star formation processcan be entirely followed through the Lagrangian point ofview of the tracers.upernova Driving. IV. The Star Formation Rate of MCs 5One of the goals of the new simulation is to describethe SN driving with a precise knowledge of the time andlocation of each SN, by resolving the formation of allmassive stars. To reduce the role of the ‘artificial’ ran-domly generated SNe, we gradually decrease their rateas the rate of the realistically generated SNe increases,so we start by reducing the random SN rate by a fac-tor of two. We continue the simulation for 10.5 Myr,with a random SN rate of only 3.12 SNe Myr − , withthe increased resolution, and with the new set of tracerparticles. At t = 55 . dx = 0 . (cid:12) (though significantly lower-mass sinkparticles are also formed). The IMF is approximately apower law, with a slope only slightly steeper than theobserved Saltpeter’s value, thus the relative number ofSNe, resulting from the explosion of the massive starsformed in the simulation, has a realistic mass and timedependence (besides a realistic spatial istribution).To follow the collapse of prestellar cores, sink parti-cles are created in cells where the gas density is largerthan 10 cm − (approximately 10 times larger than thelargest density reached by the turbulence without self-gravity), if the following additional conditions are metat the cell location: i) The gravitational potential hasa local minimum value, ii) the three-dimensional veloc-ity divergence is negative, and iii) no other previouslycreated sink particle is present within an exclusion ra-dius, r excl ( r excl = 16 dx = 0 .
12 pc in this simulation).These conditions are similar to those in Federrath et al.(2010). We have verified that they avoid the creationof spurious sink particles in regions where the gas is notcollapsing. An extensive presentation of our sink particleimplementation in Ramses can be found in Haugbølle etal. (2017).We plan to continue this simulation for ∼ −
100 Myrwith self-gravity and sink particles, in order to reach afully self-consistent solution, where the ISM turbulenceand the star formation process are driven entirely by theexplosion of massive stars formed in the simulation, in-cluding the least massive ones, of approximately 7.5 M (cid:12) ,that have a lifetime of nearly 50 Myr (Schaller et al.1992) . So far, we have reached t ≈ . (cid:12) (1408 and 2888 sink particles in the upper andlower panel respecitely) have been overplotted. Youngsinks are found inside the densest filaments, while older Neglecting even longer lifetimes for core-collapse SNe from theinteraction of intermediate-mass binaries (Zapartas et al. 2017) ones have already left their parent clouds. Some denseclusters have also cleared their surrounding gas thanksto SN explosions of their most massive members.The large number of MCs generated by this simula-tion, in combination with several hundreds of SNe fromresolved massive stars and several thousands of sink par-ticles, make up an unprecedented numerical sample tostudy the interaction of SNe with their parent clouds andto investigate the amount of clustering in space and timeof SNe of different masses. These studies are deferred toa separate paper. The analysis of the global SFR underthis fully self-consistent SN driving will also be addressedelsewhere, once the simulation will be closer to the finalgoal of ∼ −
100 Myr with self-gravity and sink par-ticles. Here, it is worth mentioning that, while produc-ing realistic SFR values within MCs, this simulation hascurrently a global SFR corresponding to a gas depletiontime of the order of 1 Gyr, consistent with global galac-tic values (Bigiel et al. 2011). Figure 3 shows the globalstar formation efficiency,
SF E , versus time (solid line),where
SF E ( t ) ≡ M s ( t ) /M box (1)and M s the total mass in sink particles. The depletiontime, t dep , is defined as: t dep ≡ M box / ( dM s ( t ) /dt ) = 1 / ( dSF E ( t ) /dt ) . (2)In Figure 3, the short-dashed and long-dashed lines show SF E versus time for depletion times t dep = 1 . t dep in the ap-proximate time intervals 61-67 Myr and 67-74 Myr. Atthe same time, the clouds selected in the simulation asdescribed in the following section exhibits a much steepertime dependence of their local SF E , as shown by theshort thin lines in Figure 3, and by the shaded areas il-lustrating the average cloud depletion time of 51 Myr.This is the first time that a value of t dep characteristicof global galactic values is derived in a simulation whereboth the star formation and its feedback are resolved, foreach individual massive stars, instead of being imposedwith subgrid-scale models. MOLECULAR CLOUD PARAMETERS
The simulation does not model the formation of ISMmolecules, so we define the clouds simply as connectedregions above a threshold density, n H , min , and refer tothem as MCs. To keep track of the effect of the valueof n H , min , we consider two values, n H , min = 200 cm − and n H , min = 400 cm − . The MC search is carried outat a uniform resolution of 512 cells, that is a spatialresolution of 0.49 pc, but the cloud properties are com-puted using all the tracer particles identified within eachcloud. Because the tracers record all the hydrodynam-ical variables interpolated at their position, and due tothe very large number of tracers in high density regions,MC properties are derived with the hydrodynamical vari-ables sampled at the highest local spatial resolution ofthe AMR grid, up to the highest resolution of 0.0076 pcin the densest regions.MCs are selected from 10 snapshots, at equal intervalsof 1.5 Myr, with the first snapshot at 4.0 Myr after theinclusion of self-gravity in the simulation. Only cloudssatisfying the following three conditions are retained: 1) Padoan et al.the cloud mass is M cl > (cid:12) , 2) the rms velocityis σ v < ), 3) the cloud does not disperse dur-ing the next 1.5 Myr, meaning that it is not doubling itseffective size in that time interval (to avoid MCs whoseproperties are evolving too rapidly). With these condi-tions, the number of MCs is reduced from 391 to 313,203 clouds with n H , min = 200 cm − and 110 clouds with n H , min = 400 cm − . The conditions allow for a bettercomparison with the SFR model that does not accountfor transient processes like SN feedback or cloud disper-sal. It may yield a sample more suitable for the compar-ison with the observations as well, because MCs stronglyaffected by SNe, or in the process of being dispersed, mayalso suffer from a strong feedback by HII regions, whichis not modeled in the simulation.In the following, we will compute for each cloud thethree non-dimensional parameters, α vir , M and β , ex-pressing the ratios of turbulent, gravitational, thermaland magnetic energies. The virial parameter, α vir , esti-mates the ratio of turbulent and gravitational energiesin a spherical cloud of radius R cl , mass M cl and one-dimensional velocity dispersion σ v (Bertoldi & McKee1992): α vir ≡ σ R cl GM cl = 403 π (cid:18) t ff t dyn (cid:19) ∼ E k E g , (3)where the dynamical time is defined as: t dyn ≡ R cl /σ v , . (4)The last equality in (3) is exact in the case of an idealizedspherical cloud of uniform density. For more realisticcloud mass distributions, the virial parameter is only anapproximation of the ratio of kinetic and gravitationalenergies. The rms Mach number is the ratio of the three-dimensional rms velocity and the sound speed, c s : M ≡ σ v , /c s , (5)and β is the ratio of gas to magnetic pressure: β ≡ P g /P m = 2 γ − ( M A / M ) , (6)where M A is the rms Alfv´enic Mach number, M A = σ v , /v A , with v A the Alfv´en velocity, γ is the adiabaticindex, and we have used the adiabatic sound speed, c s = (cid:112) γP g /ρ .To estimate the non-dimensional parameters for MCs,we need to compute the cloud mass, M cl , radius, R cl ,velocity dispersion, σ v , sound speed, c s , magnetic pres-sure, P m and thermal pressure P g . The cloud radius isdefined as the equivalent radius of the circle with areaequal to the cloud projected (along a randomly chosenaxis direction) area, A cl , R e ≡ (cid:112) A cl /π. (7)The cloud rms velocity is defined as the density-weightedone-dimensional rms velocity, and is computed from the Nearby SNe may disperse the MC (that is addressed by thethird condition for cloud selection), or only accelerate a small por-tion of the cloud, causing a strong but temporary increase of thecloud rms velocity, with very little influence on its SFR. velocity of the tracer particles: σ v ≡ (cid:34) N (cid:88) i =1 N (cid:88) n =1 ( u i,n − ¯ u i ) (cid:35) / , (8)where ¯ u i ≡ (cid:80) Nn =1 u i,n /N are the components of themean tracer particle velocity, and N is the total num-ber of tracer particles in the cloud. We choose these def-initions of radius and velocity dispersion because theyrelate directly to the observable ones, so they yield esti-mates of the virial parameter that should be comparableto the observational values. Furthermore, we have pre-viously found that the virial parameter based on thesedefinitions of radius and velocity dispersion is very closeto the actual energy ratio in realistic MCs, selected withour previous SN-driven simulation (see Sections 7 and10.5 and equation (26) in Paper I). Thus, in this workwe use the effective virial parameter, α vir , e ≡ σ R e / ( GM cl ) , (9)and the effective dynamical time, t dyn , e ≡ R e /σ v , . (10)It is important to realize that R e defined with theprojected area is the same as that defined with thecloud volume, V cl , R e , V ≡ ( V cl / (4 π/ / , only in thecase of a spherical cloud. In general, R e > R e , V ,so one should avoid estimating the mean cloud den-sity as M cl / (4 πR / SF R ff . In our cloud samples we find that (cid:104) R e /R e , V (cid:105) ≈ . t ff /t dyn , e and α vir , e becomes: t ff t dyn , e ≈ π (cid:114) α / , e ≈ . α / , e , (11)which differs by a factor of two from the coefficient de-rived from equation (3). Although we don’t know thecorrection factor in the case of real MCs, it is likelythe free-fall time (hence SF R ff ) is slightly overesti-mated as the cloud mean density is usually computedas M cl / (4 πR / c s is derived from the mass-weighted average tem-perature (the average value of the temperature associ-ated to the tracers) and σ v , = σ v √
3, using equation(8), which should yield values comparable to observa-tional estimates of MC Mach numbers. We compute β as the average of the local value of β given by the ratio ofthermal and magnetic pressure associated to each tracerparticles, with a weight m i /ρ i for each tracer, where m i and ρ i are the mass of the tracer and the local gas den-sity at the tracer position. The weight is proportional tothe volume of the gas element represented by the tracer,so we obtain a volume-averaged β . Observed values of β may be estimated from the dispersion of the polarizationangle in MCs. Due to the drop in polarization fractionwith increasing density (e.g. Padoan et al. 2001; Pelkonenet al. 2007; Ade et al. 2015), the derived β is probablycloser to a volume average value than a mass-weightedone.upernova Driving. IV. The Star Formation Rate of MCs 7 Fig. 4.—
Properties of MCs selected from the simulation with density thresholds n H , min = 200 cm − (empty circles) and n H , min = 400cm − (filled circles). The solid and dashed lines are the power-law fits to the mean values of M cl and σ v averaged within logarithmicintervals of R e (upper panels), and the values of α vir , e and β averaged inside logarithmic intervals of M cl (lower left panel) and σ v (lowerright panel). The MC properties of both samples are plotted in Fig-ure 4, where empty circles represent the clouds with n H , min = 200 cm − and filled circles the denser oneswith n H , min = 400 cm − . The two upper panels showthe mass-size and velocity-size relations. Our MCs haveequivalent radii in the range 2.6-23.5 pc, and columndensities in the range 10-80 M (cid:12) pc − . Despite the differ-ence in the average column density and size of the twocloud samples, the slopes of the relations are essentiallythe same, and are also consistent with the observationsand with results from our previous SN-driven simulation(see Papers I and II). The two lower panels of Figure 4show α vir versus M cl (left) and β versus σ v (right). Thevirial parameter spans a wide range of values, betweenapproximately 0.5 and 25, with the maximum value (andthe scatter) decreasing with increasing M cl . The gas tomagnetic pressure is in the range 0.09-1.3, with an aver-age value (cid:104) β (cid:105) = 0 .
37, and decreases with increasing σ v ,showing that the magnetic field in the clouds is ampli-fied by the cloud turbulence, as the thermal pressure doesnot vary much from cloud to cloud. However, the rmsAlfv´enic Mach number, defined as M a ≡ (cid:112) M β/ ≈ .
1, so the MC turbulence issuper-Alfv´enic (Padoan & Nordlund 1999; Padoan et al.2004; Lunttila et al. 2008, 2009).Another important non-dimensional MC parameter is the ratio, χ , of the power in compressive and solenoidalmodes of the velocity field, v , χ ≡ (cid:104) v (cid:105) / (cid:104) v (cid:105) , (12)where v c and v s are the compressive and solenoidal com-ponents of the velocity field. The two velocity compo-nents are derived with the standard Helmholtz decom-position in Fourier space, within a bounding box con-taining each MC, following the procedure in Paper II.We then compute the turbulent compressive ratio, χ t ,by subtracting the mean cloud rotation, V r , and expan-sion, V e , velocities, χ t ≡ [ (cid:104) v (cid:105) − V ] / [ (cid:104) v (cid:105) − V ] , (13)where V r and V e are derived from the mean vorticityand the mean divergence in the cloud bounding box (seeequations (1) and (2) in Paper II). As in the cloud sam-ples from our previous SN-driven simulation, we find log-normal distributions of χ t , with both mean and standarddeviations very close to our previous values, 0 . ± . n H , min = 200 cm − , and 0 . ± . n H , min = 400 cm − . The non-dimensional MC parameters will be used in Section 5to test the predictions of the turbulent fragmentationmodel. Padoan et al. Fig. 5.—
SFR per free-fall time versus effective virial parameterfor all the clouds selected from the simulation. The arrows showthe values of α vir , e of the MCs with SF R ff = 0. The error barsconnected by the solid line are the values of SF R ff averaged withinlogarithmic intervals of α vir , e . The dashed line is SF R ff ,α , the an-alytical fit to our revised PN11 model computed with the physicalparameters of the MCs from the simulation, given by equation (28).The dotted line is a function like α vir ,α , but with a smaller expo-nential coefficient, to trace roughly an upper envelope of the plot.The filled circles indicate the five MCs with M cl > × M (cid:12) (upper panel) and the six MCs with M cl > M (cid:12) (lower panel). STAR FORMATION RATE IN MOLECULAR CLOUDS
Each MC selected from the simulation is followed for1.68 Myr after the time it is identified. This time islong enough to evaluate both the time-averaged SFR andthe time variations of the SFR of individual clouds, andshort enough to avoid complications related to the iden-tification of clouds as distinct objects in the turbulentflow. MCs are continuously being formed and dispersedby the SN-driven turbulence, and they are part of aninterconnected filamentary structure that extends up tothe outer scale of the turbulence, approximately 70-100pc (Joung et al. 2009; de Avillez & Breitschwerdt 2007;Padoan et al. 2016b), so they can hardly be considered aswell-defined and long-lived isolated entities. This com-plex ever-changing nature of MCs is not a major concernfor this work, as our goal is to correlate the SFR withthe physical parameter of MCs at a given time, ratherthan to follow their long-term evolution.To define the SFR over the time interval ∆ t = 1 . SF R ≡ M tr ( t j ) − M tr ( t j + ∆ t )∆ t , (14) where, M tr ( t ) is the total mass of tracers in the cloud thathave not been accreted onto sink particles at the time t ,and t j is the time corresponding to the j -th snapshotin which the cloud is identified. The SFR per free-falltime (Krumholz & McKee 2005) is then defined as thefollowing non-dimensional SFR: SF R ff ≡ SF R / [ M tr ( t j ) /t ff ] , (15)where t ff is the free-fall time at the mean density of thecloud. To evaluate the fluctuations of SF R ff , we alsomeasure it within 14 shorter time intervals of length δt =∆ t/
14 = 0 .
12 Myr.Figure 5 shows
SF R ff versus α vir , e for our two cloudsamples (empty circles). There is a clear trend of de-creasing SFR with increasing virial parameter, thoughwith a very large scatter. The scatter increases with in-creasing α vir , e and must have primarily a random origindue to variations from cloud to cloud at fixed virial pa-rameter and time variations within each cloud, becausewe do not detect any strong dependence of SF R ff on M or β at constant α vir , e (and the scatter is much smallerin the case of the SFR predicted by our model using thecloud parameters, as shown in Figure 11). The origin ofthis scatter is discussed in Section 5.4, where we presentour model predictions. To illustrate the amount of timevariations, we show in Figure 6 the maximum and mini-mum values of SF R ff measured in the 14 time intervals δt = 0 .
12 Myr in each cloud. The characteristic scat-ter of
SF R ff due to time variations is over one order ofmagnitude.The mean values of SF R ff in logarithmic bins of α vir , e are shown by the error-bar symbols in Figure 5, wherethe size of the error bars is equal to the standard er-ror of the mean (the rms divided by square root of thethe number of data points in the bin). Despite the siz-able scatter, these mean values follow very closely theanalytical fit to our model (see Section 5.3) applied toall the MCs selected from the simulation, SF R ff ,α =0 . − . α / , e ). This agreement is nearly equallygood for both cloud samples, with the lower density one( n H , min = 200 cm − ) yielding only slightly smaller meanvalues (though always within approximately one stan-dard error of the mean from the analytical fit of themodel). The agreement with the model breaks downonly at large values of α vir , e , where the statistical signifi-cance of the mean values of SF R ff is low, due to the verysmall number of data points per bin (the standard errorof the mean can still be small, if most points are close toeach other by chance).The arrows in Figure 5 indicate the values of α vir , e ofthe clouds where SF R ff = 0. The upper panel showsthat a significant fraction of clouds in the sample with n H , min = 200 cm − (approximately 18%) have SF R ff =0, even at relatively low values of α vir , e . The fractiondrops to 5% in the case of n H , min = 400 cm − (lowerpanel). We find that only relatively small clouds have SF R ff = 0. Figure 7 shows that SF R ff = 0 only inclouds with M cl < × M (cid:12) or M cl < × M (cid:12) inthe lower and higher n H , min samples respectively. Thelower envelope of the plots in Figure 7 rises sharply withincreasing cloud mass, so that the scatter in SF R ff israther small for the most massive clouds. The averagevalues for the 5 and 6 most massive clouds in the twoupernova Driving. IV. The Star Formation Rate of MCs 9 Fig. 6.—
Time variation of
SF R ff in each cloud, with cloudsordered by increasing value of their time-averaged SF R ff (solidblack line). Each error bar extends between the maximum andminimum SFR of a cloud, among the values measured within 14time intervals of ∆ t/
14 = 0 .
12 Myr size. Values of 10 − indicatethat the minimum SFR is zero. Clouds that never form stars duringthe whole interval ∆ t = 1 .
68 Myr are not shown in the figure (35clouds in the catalog with n H , min = 200 cm − and 6 clouds for n H , min = 400 cm − .) samples (shown as filled circles in Figures 5 and 7) are SF R ff = 0 . ± .
01 ( n H , min = 200 cm − ) and SF R ff =0 . ± .
02 ( n H , min = 400 cm − ).The probability distributions of the SFR are shown inFigure 8. Although they both peak at SF R ff ≈ . SF R ff = 0 . ± .
03 for n H , min = 200 cm − and SF R ff = 0 . ± .
04 for n H , min = 400 cm − . Thesevalues are comparable to observational estimates of theSFR in MCs, as discussed below in Section 6.In Figure 9, we plot the mass fraction turned into starsin the time interval ∆ t = 1 .
68 Myr after the cloud is iden-tified, which is given by
SF R × ∆ t , versus the initial starformation efficiency at the time the cloud is identified, SF E , given by SF E ≡ M tr , s ( t j ) / ( M tr , s ( t j ) + M tr ( t j )) , (16)where M tr , s ( t j ) is the total mass of the tracer particlesinside the cloud and already accreted onto sinks at thetime t j when the cloud is identified, and M tr ( t j ) is thetotal mass in tracers found in the same cloud and still inthe gas phase.We can express the initial SFE as SF E = SF R × t cl ,where SF R is the average SFR up to the time the cloudis identified ( t < t j ), and t cl is the cloud age. Althoughonly some of the MCs have an approximately constantSFR during the time interval ∆ t , if we assumed thatthe average SFR prior to the cloud identification were Fig. 7.—
SFR per free-fall versus cloud mass for the same MCsfrom the simulation as in the previous figures.
Fig. 8.—
Probability distributions of
SF R ff of our two samplesof MCs extracted from the simulation. equal to the average SFR after the cloud identification, SF R = SF R , we could derive the age of a cloud, t cl ,at the time when it is selected, as t cl ≡ SF E /SF R .Constant values of t cl are shown in Figure 9 by dashedlines. The upper envelope of the plot in Figure 9 is verywell described by a line of constant cloud age, t cl ≈ . SF R = 0and
SF E > Fig. 9.—
Cloud mass fraction converted into stars in the timeinterval during which we follow the MCs, ∆ t = 1 .
68 Myr, versusthe cloud initial SFE,
SF E . Assuming that the SFR prior tothe time of the cloud selection, SF R , was constant and equalto the average one measured during ∆ t after the cloud selection, SF R = SF R , we derive an estimate of the MC age at the timeof selection, t cl ≡ SF E /SF R , shown by the dashed lines for agesof 1.5, 6.0 and 24 Myr. As in previous figures, the filled circles arethe five MCs with M cl > × M (cid:12) (upper panel) and the sixMCs with M cl > M (cid:12) (lower panel). The arrows indicate thevalues of SF E for clouds with SF R = 0. of constant SFR does not apply and their SFR has de-clined with time,
SF R > SF R , rather than their agebeing much larger than 20 Myr.The MC lifetime was studied in Paper I, using ourprevious SN-driven simulation, where we found that itwas approximately four dynamical times on average, t life /t dyn = 4 . ± .
7, resulting in cloud lifetimes con-sistent with those estimated for MCs in the Large Mag-ellanic Cloud (Kawamura et al. 2009). Using the aver-age cloud dynamical times in our two cloud samples, wewould estimate an average cloud lifetime t life = 14 . ± . n H , min = 200 cm − and t life = 11 . ± . n H , min = 200 cm − . Most clouds in the two sampleshave smaller estimated ages than such lifetimes, suggest-ing that most of the ages derived from the assumption ofconstant SFR are reasonable. THE TURBULENT FRAGMENTATION MODEL
The SFR can be modeled as the result of the gas frag-mentation by the supersonic turbulence (Padoan 1995;Krumholz & McKee 2005; Padoan & Nordlund 2011;Federrath & Klessen 2012; Hopkins 2013). Given theprobability density function (PDF) of the gas density insupersonic MHD turbulence and a critical density for col-lapse, the integral of the PDF above the critical density,divided by an appropriate time-scale, provides an esti-mate of the SFR. Different models differ from each otherin modeling the PDF, in the definition of the critical den- sity and in the choice of the time-scale, and were recentlyreviewed in Padoan et al. (2014). In this work, we onlyconsider the model by Padoan & Nordlund (2011, PN11hereafter), which we slightly revise, motivated both byphysics considerations and by the numerical results. Inthe following subsections, we first summarize the PN11model and present our revision; we then apply the re-vised model to the clouds selected from the simulationand discuss the origin of the scatter in the
SF R ff – α vir relation. The PN11 Star Formation Rate Model
The model depends on the three non-dimensional pa-rameters, α vir , M and β , that express the ratios of tur-bulent kinetic, gravitational, thermal and magnetic en-ergies. These parameters were defined and computed inSection 3 for all the clouds of our two samples. In themodel revision described below, we introduce a fourthnon-dimensional parameter, χ t , also defined in Section3. The critical density is defined as the external density ofa critical Bonnor-Ebert sphere of size equal to the char-acteristic thickness of the postshock gas in the turbulentflow. Using the MHD jump conditions of an isothermalshock, PN11 derive the following expression: ρ cr ρ = 0 . θ − α vir M (1 + 0 . β − ) (1 + β − ) , (17)where ρ is the mean gas density, and θ is the fractionof the cloud diameter corresponding to the characteristicsize of the largest compressive motions in the turbulentflow, with θ = 0 .
35 in PN11.The density PDF is assumed to be a lognormal distri-bution, p ( x ) dx = 1 x (2 πσ ) / exp (cid:20) − (ln x + σ / σ (cid:21) dx, (18)where x = ρ/ρ is the gas density normalized to the meanand σ is the standard deviation of ln( x ) that depends onboth M and β , σ ≈ ln (cid:104) b M ) (1 + β − ) − (cid:105) , (19)corresponding to the following standard deviation for x , σ x ≈ (1 + β − ) − / b M . (20)In PN11, b = 0 . β = 0 .
39, based on numerical es-timates and fitting of the PDF. In Section 3, we haveshown that the average value of beta for the clouds se-lected from the simulation is (cid:104) β (cid:105) = 0 .
39, while the aver-age value of b , using its dependence on χ t given belowand proposed in Paper II, and the distribution of χ t forour clouds, is (cid:104) b (cid:105) = 0 .
48. Thus, the values of b and β adopted in PN11 were very close to the ones derived inthis work for a realistic sample of MCs.Assuming that a fraction (cid:15) of the mass fraction abovethe critical density is turned into stars in a free-falltime of the critical density, t ff , cr = (3 π/ (32 Gρ cr )) / , thestar formation rate per free-fall time (the mass fractionupernova Driving. IV. The Star Formation Rate of MCs 11 Fig. 10.—
SFR per free-fall time versus virial parameter pre-dicted by our revision of the PN11 model, for five combinationsof four values of M and three values of β , assuming (cid:15) = 0 . b = 0 .
48 (the latter is the average value from our cloud samples).The dashed line shows the analytical fit to our model applied to theMCs from the simulation,
SF R ff ,α (see Section 5.3 and Figure 11),and the dotted line the same function, but with a smaller exponen-tial coefficient, that was used in Figure 5 to trace an approximateupper limit of the SF R ff - α vir , e relation. turned into stars in a free-fall time) is given by: SF R ff , PN11 = (cid:15) t ff t ff , cr (cid:90) ∞ x cr x p ( x ) dx = (cid:15) x / (cid:18) (cid:20) σ − x cr )2 / σ (cid:21)(cid:19) (21)where t ff = (3 π/ (32 Gρ )) / is the free-fall time of themean density and x cr = ρ cr /ρ . Revision of the PN11 SFR Model
The first modification to our model is the choice ofthe time-scale that defines the SFR. While in PN11 weassumed that the timescale was the free-fall time of thecritical density, t ff , cr , here we choose the time of forma-tion of the high-density tail of the PDF, t PDF , which wedefine as the lifetime of compressions responsible for thecharacteristic postshock density used in our derivation ofthe critical density, and we also assume R cl = R e , t PDF ≡ θ R e /σ v , = 2 θ t dyn , e . (22)Thus, the coefficient of SF R ff is t ff /t PDF instead of t ff /t ff , cr . The free-fall time of the critical density is gener-ally too short to assume that the high-density tail of thePDF is maintained despite the collapse of the dense gas.In our cloud samples, (cid:104) t PDF /t ff , cr (cid:105) ≈ .
8, so the averageSFR is set by t PDF rather than t ff , cr , while t ff , cr (and theunrevised PN11 model) sets the maximum value of theSFR, when the PDF tail is fully sampled.The ratio t ff /t PDF can be related to α vir , e using equa-tion (11): t ff t PDF ≈ π θ (cid:114) α / , e ≈ . θ − α / , e , (23)so our revised expression for the SFR per free-fall time is: SF R ff , MHD = 0 . (cid:15) θ − α / , e (cid:18) (cid:20) σ − x cr )2 / σ (cid:21)(cid:19) . (24)The second modification is in the expression for thecritical density. Because it is based on the characteristicpostshock density, it makes sense to consider the prod-uct b M instead of M /
2, in analogy with the derivationof the standard deviation of the density PDF. Further-more, because we can measure the ratio, χ t , of the powerin compressive and solenoidal modes for each cloud, werelate the parameter b to χ t by assuming, as in PaperII, that b M is the rms Mach number of the compressivepart of the velocity field, which gives b = (cid:112) χ t / (1 + χ t ) . (25)With this final modification, b and M only appear intheir product, so a change in b is equivalent to a changein M .The revised version of the critical density is: x cr = 0 . θ − α vir , e χ t (1 + χ t ) M (1 + 0 . β − ) (1 + β − ) , (26)This revised model is illustrated in Figure 10, for fivecombinations of four values of M and three values of β , assuming (cid:15) = 0 . b = 0 .
48, the average value ofour cloud sample. The dashed line shows the analyticalfit to the model applied to the MCs from the simula-tion,
SF R ff ,α = 0 . − . α / , e ) (see Section 5.3),as in Figure 5. The reference model (medium-size filledcircles) with the parameter values corresponding to thepeak (not the average) of the probability distributionsof M and β , M = 7 and β = 0 .
2, is nearly indistin-guishable from the analytical fit,
SF R ff ,α . The figurealso shows that a variation of M by a factor f corre-sponds approximately to a variation of β by a factor f .This can be easily seen in the limit of β (cid:28)
1, whereboth coefficients containing β in equations (17) and (19)become ≈ β , so β and M only appear in the product β M , which explains the observed dependence on theseparameters in Figure 10. Thus, in the limit of β (cid:28) α vir , e , and the effective Mach number, M e , M e ≡ b M β / = (2 /γ ) / b M A (27)The second equality, derived from the definition of β ineq. (6), shows that, in the limit β (cid:28)
1, the effectiveMach number is the compressive part of the Alfv´enicMach number. Furthermore, at a value of α vir , e ≈ − SF R ff , MHD is nearly independent of M e , as shownin Figure 10, while it increases with increasing M e for α vir , e (cid:38) M e for α vir , e (cid:46) Model Predictions for the Numerical MC Samples
We apply the revised PN11 model to the physicalparameters of the numerical samples of MCs derivedin Section 3. The dependence of the predicted SFR,
SF R ff , MHD , on α vir , e is shown in Figure 11 (large emptycircles). Despite the wide range of parameter values,the predicted SF R ff , MHD - α vir , e relation has a very smallscatter, compared with that of the corresponding numer-ical relation in Figure 5. This is partly due to the anti-correlation between M and β , as shown in Figure 12,resulting in a small scatter in M e . The average andstandard deviation of M e are 2 . ± . . ± . n H , min samples respectively, whichexplains the small scatter in Figure 11.Because the predicted SF R ff , MHD - α vir , e relation doesnot show a significant dependence on n H , min , and givenits small scatter, we provide an analytical fit inspired byour earlier results from a large set of idealized turbulencesimulations (Padoan et al. 2012): SF R ff ,α = 0 . − . α / , e ) , (28)shown by the dashed line in Figure 11. In Section 4 weshowed that SF R ff ,α is also an excellent fit to the val-ues of SF R ff from the simulation averaged within log-arithmic intervals of α vir , e (Figure 5). Thus, the modelpredicts successfully the average SFR, but without de-scribing its scatter. We also showed, in Section 5.2,that SF R ff ,α is nearly indistinguishable from the re-vised model, SF R ff , MHD , with a fixed set of parameters, (cid:15) = 0 . b = 0 . β = 0 . M = 7 (see Figure 10).As shown in Figure 12, M and β are correlated with α vir , e , in such a way that the product M β / is nearlyindependent of α vir , e . Furthermore, we find that b is alsoindependent of α vir , e , so the effective Mach number, M e ,is also approximately independent of α vir , e , and varia-tions of M e from cloud to cloud cannot modify the rela-tion SF R ff , MHD – α vir , e , but only contribute to its (small)scatter. This is the reason why a single set of parametervalues provides a relatively good fit to the average SFRof all clouds.In Padoan et al. (2012), using a large set of AMR sim-ulations, we derived a SFR law that depended only onthe ratio t ff /t dyn or, equivalently, on α vir , SF R ff , P12 = (cid:15) exp( − . t ff /t dyn ). Using equation (3) and assumingthat the efficiency factor is (cid:15) = 0 .
5, we would ob-tain
SF R ff , P12 = 0 . − . α / ), so the coefficientswould be different from those of the analytical fit in thiswork, SF R ff ,α . However, the value of t ff /t dyn in Padoanet al. (2012) was the average over the whole computa-tional volume, while most of the star formation in thosesimulations occurred in dense clumps where the localvalue of t ff /t dyn (or α vir ) may have been smaller than themean value. The results in Padoan et al. (2012) are notnecessarily inconsistent with the current results , but re- They would be consistent if, for example, the star-formingclumps in the turbulence simulations had a virial parameter onaverage approximately 30% smaller than the global one.
Fig. 11.—
SFR per free-fall time versus effective virial parameterpredicted by our model (equation (24)) for the physical parametersof the MCs extracted from the simulation (large empty circles).The filled circles show the five (upper panel) and six (lower panel)most massive MCs, as in previous figures. The dashed line is thefit to the clouds with α vir , e < SF R ff ,α , given in equation (28).The dotted line is the same function, but with a smaller exponentialcoefficient, that was used in Figure 5 to trace an approximate upperlimit of the SF R ff - α vir , e relation. The small empty circles areupper limits for each cloud predicted by the model (see text inSection 5.3). lating the SFR in those idealized turbulence simulationsto the SFR in real MCs is an open problem. The need fora more realistic numerical sample of star-forming cloudswas one of the motivations for the current work. On theother hand, results from a general cubic region of theISM, as in Padoan et al. (2012), are more relevant if thegoal is to develop a subgrid model for star formation forgalaxy formation simulations (e.g. Semenov et al. 2016). The Scatter in the
SF R ff – α vir , e Relation
The comparison of Figures 5 and 11 shows that thescatter in the values of
SF R ff at constant α vir , e fromthe simulation is much larger than the predicted scatterfrom our revised PN11 model, so the variations in b , β or M from cloud to cloud can explain only a small fractionof the scatter. However, large variations from cloud tocloud and time variations in individual clouds are notinconsistent with the physical assumptions of the model.The average ratio between the formation time of densestructures in the revised model and their collapse timeis (cid:104) t PDF /t ff , cr (cid:105) ≈ .
8, where the value is computed overboth cloud samples. Thus, there could be periods ofhigh SFR, lasting for a time of order (cid:104) t ff , cr (cid:105) ≈ . (cid:104) t PDF (cid:105) ≈
Fig. 12.—
Ratio of gas pressure to magnetic pressure (upperpanel) and rms Mach number (lower panel) versus effective virialparameter for the MCs selected in the simulation. The solid anddashed lines are the power-law fits to the mean values of β and M averaged inside logarithmic intervals of α vir , e , for the cloud sampleswith n H , min = 200 (solid lines) and 400 (dashed lines) cm − . level. This intermittent behavior is more likely to occurin regions where the SFR is very low, such as in smallclouds and for high values of α vir , e (which are also morelikely to occur in lower mass clouds), because of the smallnumber of collapsing cores. In very large clouds and atsmall values of α vir , e , the number of cores is very largeand time variations of the high density tail of the PDFshould have a smaller amplitude. In Figure 5, the scatterindeed decreases towards smaller values of α vir , e .The model is consistent with finding clouds with nostar formation, particularly at large α vir , e , if the PDFtail is completely depleted. The maximum SFR can beestimated from the integral of the PDF above the criti-cal density, assuming the PDF tail is fully sampled, di-vided by the collapse time, so the coefficient of SF R ff is t ff /t ff , cr , as in PN11, instead of t ff /t PDF as in the revisedmodel. The estimated values of the maximum SFR for allthe clouds in our samples are shown in Figure 11 (smallempty circles). They are almost a factor of three to fourlarger than the predicted mean values (very close to theprediction of the PN11 model) and follow approximatelythe dotted line that was shown to be an approximateupper envelope of the
SF R ff - α vir , e relation in Figure 5.At α vir , e <
1, the maximum values are systematicallyabove the dotted line, while the
SF R ff values from thesimulation are all below the line, consistent with our ex-pectation that the scatter caused by time variations ofthe high-density tail of the PDF should be smaller forsmall α vir , e .To evaluate the amplitude of the deviations in thehigh-density tail of the cloud density PDF relative to the Fig. 13.—
Measured versus predicted mass fraction above thecritical density for every cloud in the two numerical catalogs. M cr is the cloud mass above the critical density measured in the sim-ulation at the time of cloud selection; M cr , MHD is the cloud massabove the critical density predicted by the model.
Fig. 14.—
SFR per free-fall time versus effective virial param-eter for a subset of the clouds selected from the simulation (fromboth cloud catalogs) harboring more than 100 sink particles (filledcircles). The model prediction,
SF R ff , MHD , for each cloud is alsoshown (empty squared symbols). For reference, we also show thesame dashed and dotted curves as in Figures 5 and 11. The scatterin
SF R ff for this cloud subsample is significantly reduced relativethat of the full samples shown in Figures 5, and the deviations fromthe model predictions are on average less than 50%. model prediction, we compute, for each cloud, the massfraction above the predicted critical density based on thecloud density field from the simulation, M cr /M cl , where M cr is the mass above the critical density. We then com-pare this mass fraction with that based on the densityPDF predicted by the model for the physical parametersof the cloud, M cr , MHD /M cl , where M cr , MHD is the mass4 Padoan et al.above the critical density in the theoretical PDF. Thecomparison is shown in Figure 13 for both cloud catalogs.Because larger values of these mass fractions correspondto higher values of the
SF R ff , Figure 13 exhibits thesame trend of increasing scatter with decreasing SF R ff as in Figure 5. The scatter is relatively large, though notenough to explain the corresponding scatter in SF R ff inFigure 5. Furthermore, many clouds still deviate sig-nificantly from the predicted SFR even if we apply themodel with the mass fraction above the critical densitydirectly measured in the clouds, M cr /M cl , rather thanthe one predicted by the model, M cr , MHD /M cl . Thus,variations of the density PDF cannot be the only expla-nation for the scatter in the relation between SF R ff and α vir , e found in the simulation.Random variations in SF R ff due to the finite num-ber of collapsing cores in each cloud certainly contributeto the scatter, particularly in clouds with high values of α vir , e that contain a relatively small number of sink par-ticles. This naturally explains the increase in the scat-ter with increasing α vir , e in Figure 5. We can illustratethe role of low-number statistics by selecting a subsam-ple of clouds containing a large number of sink particles, N sink , which should significantly reduce the effect of ran-dom fluctuations of the specific physical realization of theturbulence in each cloud. In Figure 14 we show the caseof clouds with N sink > SF R ff , MHD , for each cloud (empty squared sym-bols). The scatter is very much reduced compared withFigure 5, and the values of
SF R ff from the simulationdeviate from the model predictions by less than 50% onaverage. In Figure 14, the scatter does not show theclear dependence on α vir , e seen in Figure 5, as we haveexcluded the clouds with low sink-particle numbers. Theeffect of random variations in SF R ff may be somewhatlarger in our simulation than in nature, due to the in-completeness of our stellar IMF below approximately 5-10 M (cid:12) . This will be tested in future works using zoom-insimulations of individual clouds, where the increase res-olution will yield a larger number of collapsing cores anda complete IMF down to a fraction of a solar mass.Additional scatter, contributing to the total one in Fig-ure 5, must arise from the lack of statistical equilibriumin the clouds, due to the specific boundary conditions ofeach cloud and the limited cloud lifetime. The averagecloud lifetime is rather short, t life ≈ . t dyn (see PaperI), and MCs are hardly isolated objects: they are partof a complex filamentary network where mass accretiononto the clouds may not be negligible, and the cloudstructure is ever changing, with the cloud ultimately be-ing dispersed. While the MC turbulence is driven by SNewith an effective outer scale of order 70-100 pc (see PaperI), nearby SNe can also affect the clouds directly, tem-porarily changing the relation between SF R ff and α vir , e (for example causing a sudden increase in the cloud rmsMach number, while only affecting a relatively small frac-tion of the cloud mass and thus not modifying the SFRsignificantly.)As explained in Section 4, we computed the SFR asan average over a time ∆ t = 1 .
68 Myr. Although thescatter in the SFR may be reduced by averaging overa significantly longer time (the average dynamical timeof the clouds in our samples is nearly 3.0 Myr, so the average lifetime may be of order 13 Myr), we did not tryto define the SFR over larger ∆ t . This is because largevariations from cloud to cloud are also found from theobservations, thus it is important to avoid measuring theSFR in the simulation in a way that cannot be relatedto the observational SFR estimates. As discussed below,the SFR is measured in nearby MCs by counting thenumber of protostars, and assuming a typical protostellarlifetime of 2 Myr, which is probably a reasonable estimateof the duration of class II protostars (Evans et al. 2009;Spezzi et al. 2008). Even when the SFR is measuredin a sample of clusters selected by free-free emission ofionized gas, the typical cluster lifetime is probably of theorder of 2 Myr, at least if the star formation process isvery rapid, because the free-free emission is expected tostrongly decay in a time of order 4 Myr, essentially thelifetime of the massive stars that contribute the most tothe ionized flux (Murray 2011). Thus, our value of ∆ t iscomparable to the characteristic age of the young stellarpopulations used in the determination of the SFR in realMCs. Furthermore, as explained in Section 4, we avoidlonger ∆ t also to limit uncertainties related to the cloudidentification, such as the possibility of a significant massaccretion (which we neglect) or the chance that a cloudis completely dispersed during the time we measure theSFR.Given this unavoidable scatter in SF R ff at constant α vir , e , we conclude that SFR models should be comparedwith simulations or observations by computing ensem-ble averages with large cloud samples. This has beenachieved with our simulation, allowing us to show that,despite the scatter, the values of SF R ff averaged in loga-rithmic bins of α vir , e closely follow the revised model (seeFigure 5). Observational efforts should also aim at com-piling very large samples of MCs with estimated valuesof SF R ff and cloud parameters. COMPARISON WITH OBSERVATIONS
The average SFR in MCs can be deduced from theglobal SFR in the Galaxy. For example, Krumholz &Tan (2007) derived a value of
SF R ff = 0 .
02, and Mur-ray (2011) a value approximately three times smaller,
SF R ff = 0 . SF R ff = 0 . − .
24, and, more importantly, a verywide range of values, from 0.001 to 0.59. Vutisalchavakulet al. (2016) selected star-forming regions based on acatalog of HII regions (Anderson et al. 2014), but de-rived much lower SFRs, based on the 22 µ m band ofthe WISE satellite (Wright et al. 2010), with an averagevalue of SF R ff = 0 . Fig. 15.—
SFR per free-fall time versus effective virial param-eter for the MC sample by Evans et al. (2014) (filled circles) andfor the MCs selected from our simulation with density threshold n H , min = 400 cm − (empty circles). The dashed line is the sameanalytical fit to our model prediction, SF R ff ,α , as in previous fig-ures. The dotted line is our model prediction for a single set ofparameters, using the average Mach number value in the observa-tional sample, M = 9 . b and β as the average values from the simulation, b = 0 . β = 0 .
37 (because they are unknown for the observed clouds).The green shaded area shows the estimated
SF R ff and α vir , e forthe CMZ clouds, from Federrath et al. (2016) and Barnes et al.(2017). velocity dispersion may be significantly affected by thefeedback from massive stars, especially in the case of themost active star-forming clouds. This would produce adependence of the virial parameter on the SFR, ratherthan probing the dependence of the SFR on the virialparameter.To avoid the large uncertainties in the properties ofdistant clouds associated with HII regions, we compareour numerical results with the observations using onlydata from well-studied nearby MCs. Evans et al. (2014)derived MC properties and SFRs for 29 clouds from thec2d (Evans et al. 2003, 2009) and Gould Belt (Dunhamet al. 2013) Spitzer legacy programs. They found anaverage value of
SF R ff = 0 .
018 (0.016 including the fourclouds with
SF R ff = 0), by counting all the protostarsin the clouds, and assuming a mean stellar mass of 0.5M (cid:12) and a timescale of 2 Myr for Class II protostars.In Figure 15, we plot SF R ff versus α vir , e for thesenearby MCs (filled circles), and the α vir , e values for thefour clouds with SF R ff = 0 (large arrows at α vir , e > n H , min = 400 cm − (empty circles). Including the cloudswith SF R ff = 0, SF R ff decreases with increasing α vir , e ,though not as rapidly as in the simulation. However, theaverage Mach number of these clouds is 9.4, comparedwith the value of 7.7 from our numerical samples (usingthe same cloud mean temperature in both cases), so the SF R ff - α vir , e relation is expected to be a bit shallowerfor the observed clouds than for the numerical samples.The dotted line in Figure 15 shows the prediction of ourrevised model with M = 9 .
4, while the dashed line is theanalytical fit to the model applied to the clouds from thesimulation, as in previous figures. The number of cloudsin the observational sample is too low to define a clearrelation between
SF R ff and α vir , e , but the cloud SFRs are approximately distributed around the model predic-tion shown by the dotted line. Furthermore, the valuesof SF R ff of the clouds from the simulation appear tobe consistent with those in nearby MCs. The scatter in SF R ff increases with increasing α vir , e both in the sim-ulation and in the observations. The nearby MCs with α vir , e (cid:38) SF R ff = 0 − . SF R ff and α vir , e are relatively well determined for theCMZ (Barnes et al. 2017). The SFR value derived byBarnes et al. (2017), SF R ff = 0 . − .
04, associatedwith the value α vir , e = 4 . ± . SUMMARY AND CONCLUSIONS
We have studied the SFR as a function of cloud pa-rameters by generating a large sample of realistic MCs,formed ab initio in a simulation of SN-driven turbulencein an ISM region of 250 pc size. This simulation is acontinuation of our previous SN-driven experiment pre-sented in Papers I, II and III, where we had alreadydemonstrated that we could recover MC properties con-sistent with the observations. In this work, thanks to asignificant increase of the spatial resolution and the in-troduction of sink particles, we can further test if theSFR in the clouds is consistent with the observations aswell. Although the global SFR in the simulation will beaddressed elsewhere, we have anticipated that it corre-sponds to a depletion time of order 1 Gyr, in agreementwith global galactic values, which is also an importanttest for the simulation. The main results of this workare summarized in the following.1. The SFR per free-fall time in the MCs selected fromthe simulation follows a broad probability distribu-tion, with a peak at
SF R ff ≈ . SF R ff in the simulation decreases withincreasing α vir , e .3. The SF R ff - α vir , e relation from the simulation has alarge scatter that is not explained by cloud to cloudvariations of the other non-dimensional parameters b , β and M . This scatter is most likely due to acombination of time variations in the high-density6 Padoan et al.tail of the gas density PDF, random fluctuations of SF R ff in clouds with a low number of stars (sinkparticles), and a lack of statistical equilibrium ofthe MC turbulence, due to the transient nature ofthe MCs.4. The PN11 model has been revised, with the mostimportant modification being the choice of thetimescale. While in PN11 we chose the free-falltime of the critical density, t ff , cr , in the revisedmodel we assume it is the timescale of formationof the characteristic post-shock structures responsi-ble for the high-density tail of the gas density PDF, t PDF = 2 θ t dyn , because t ff , cr < t PDF . This choiceresults in a SFR a few times lower than in the PN11model, for characteristic parameters of MCs.5. Applied to the MCs from the simulation, the re-vised model results in a
SF R ff - α vir , e relation witha rather small scatter, which is fit well by the re-lation SF R ff ,α = 0 . − . α / ). This relationis consistent with our previous result in Padoanet al. (2012), where we had already concluded that SF R ff depends primarily on the virial parameter.6. The values of SF R ff in the MCs selected fromthe simulation, averaged in logarithmic intervalsof α vir , e , follow the prediction of our revised PN11model, SF R ff ,α (at least for α vir , e (cid:46) N sink > SF R ff is expected to be small.7. The SFR measured in well studied nearby MCsfrom direct counts of protostars is consistent withthe SFR values of the clouds in the simulation andwith our model predictions. As in the simulation,the scatter is large also in the observational valuesof SF R ff , and increases towards larger values of α vir , e .8. The SFR in the CMZ, estimated with the aid ofa study of the physical properties of the “Brick” cloud and of an orbital model of the clouds in theCMZ, is also consistent with our numerical resultsand our theoretical predictions.One of the most valuable findings of this investigationis the large scatter of SF R ff , including both large varia-tions from cloud to cloud, even at comparable values of α vir , e , and time variations within individual clouds. Thisresult has been achieved thanks to the very large numberof clouds formed ab initio in the simulation, each con-tributing to an ensemble of realistic initial and bound-ary conditions with probability distributions that mayclosely match those of real MCs. The instantaneous valueof the SF R ff of a given cloud cannot be accurately pre-dicted based on the cloud physical parameters, and caneven deviate significantly from the prediction of the the-oretical model, because of the ever-changing nature ofMCs embedded in the complex filamentary structure ofthe cold ISM, and the continuous driving by SNe, withoccasional explosions in close proximity or even withinthe cloud. Despite this chaotic nature of the SFR inMCs, we have found that an ensemble average obtainedby selecting a very large number of clouds from manysnapshots of the simulation yields values of SF R ff withinlogarithmic intervals of α vir , e that nicely fit the modelprediction. We conclude that our simulation providessupport for the theoretical framework of turbulent frag-mentation, while exposing the chaotic and unpredictablenature of the star formation process.We thank the anonymous referee for several usefulcomments and corrections that helped us improve thepaper. Computing resources for this work were pro-vided by the NASA High-End Computing (HEC) Pro-gram through the NASA Advanced Supercomputing(NAS) Division at Ames Research Center. PP acknowl-edges support by the Spanish MINECO under projectAYA2014-57134-P. TH is supported by a Sapere AudeStarting Grant from The Danish Council for Indepen-dent Research. Research at Centre for Star and PlanetFormation is funded by the Danish National ResearchFoundation. REFERENCESAde, P. A. R., Aghanim, N., Alina, D., Alves, M. I. R.,Armitage-Caplan, C., Arnaud, M., Arzoumanian, D., Ashdown,M., Atrio-Barandela, F., & et al. 2015, A&A, 576, A104Anderson, L. D., Bania, T. M., Balser, D. S., Cunningham, V.,Wenger, T. V., Johnstone, B. M., & Armentrout, W. P. 2014,ApJS, 212, 1Bertoldi, F., & McKee, C. F. 1992, ApJ, 395, 140Bigiel, F., Leroy, A. K., Walter, F., Brinks, E., de Blok, W. J. G.,Kramer, C., Rix, H. W., Schruba, A., Schuster, K.-F., Usero,A., & Wiesemeyer, H. W. 2011, ApJ, 730, L13Bonnell, I. A., Dobbs, C. L., & Smith, R. J. 2013, MNRAS, 430,1790de Avillez, M. A., & Breitschwerdt, D. 2007, ApJ, 665, L35Dekel, A., Sari, R., & Ceverino, D. 2009, ApJ, 703, 785Dobbs, C. L. 2015, MNRAS, 447, 3390Dobbs, C. L., & Pringle, J. E. 2013, MNRAS, 432, 653Dunham, M. M., Arce, H. G., Allen, L. E., et al. 2013, AJ, 145, 94Evans, N. J., II, Allen, L. E., Blake, G. A., et al. 2003, PASP,115, 965Evans, N. J., II, Dunham, M. M., Jørgensen, J. K., et al. 2009,ApJS, 181, 321-350 Evans, II, N. J., Heiderman, A., & Vutisalchavakul, N. 2014, ApJ,782, 114Federrath, C., Banerjee, R., Clark, P. C., & Klessen, R. S. 2010,ApJ, 713, 269Federrath, C., & Klessen, R. S. 2012, ApJ, 761, 156Federrath, C., Rathborne, J. M., Longmore, S. N., Kruijssen,J. M. D., Bally, J., Contreras, Y., Crocker, R. M., Garay, G.,Jackson, J. M., Testi, L., & Walsh, A. J. 2016, ArXiv e-printsForbes, J. C., Krumholz, M. R., Burkert, A., & Dekel, A. 2014,MNRAS, 438, 1552Fromang, S., Hennebelle, P., & Teyssier, R. 2006, A&A, 457, 371Heiles, C., & Troland, T. H. 2005, ApJ, 624, 773Heyer, M. H., Brunt, C., Snell, R. L., Howe, J. E., Schloerb, F. P.,& Carpenter, J. M. 1998, ApJS, 115, 241Heyer, M. H., Carpenter, J. M., & Snell, R. L. 2001, ApJ, 551, 852Hopkins, P. F. 2013, MNRAS, 430, 1653Iffrig, O., & Hennebelle, P. 2015, A&A, 576, A95Joung, M. R., Mac Low, M.-M., & Bryan, G. L. 2009, ApJ, 704,137Kauffmann, J., Pillai, T., Zhang, Q., et al. 2016, arXiv:1610.03499 upernova Driving. IV. The Star Formation Rate of MCs 17upernova Driving. IV. The Star Formation Rate of MCs 17