The Global Evolution of Giant Molecular Clouds II: The Role of Accretion
Nathan J. Goldbaum, Mark R. Krumholz, Christopher D. Matzner, Christopher F. McKee
aa r X i v : . [ a s t r o - ph . GA ] J un Draft version April 18, 2018
Preprint typeset using L A TEX style emulateapj v. 11/10/09
THE GLOBAL EVOLUTION OF GIANT MOLECULAR CLOUDS II:THE ROLE OF ACCRETION
Nathan J. Goldbaum , Mark R. Krumholz , Christopher D. Matzner , Christopher F. McKee Draft version April 18, 2018
ABSTRACTWe present virial models for the global evolution of giant molecular clouds. Focusing on the presenceof an accretion flow, and accounting for the amount of mass, momentum, and energy supplied byaccretion and star formation feedback, we are able to follow the growth, evolution, and dispersal ofindividual giant molecular clouds. Our model clouds reproduce the scaling relations observed in bothgalactic and extragalactic clouds. We find that accretion and star formation contribute roughly equalamounts of turbulent kinetic energy over the lifetime of the cloud. Clouds attain virial equilibrium andgrow in such a way as to maintain roughly constant surface densities, with typical surface densities oforder 50 – 200 M ⊙ pc − , in good agreement with observations of giant molecular clouds in the MilkyWay and nearby external galaxies. We find that as clouds grow, their velocity dispersion and radiusmust also increase, implying that the linewidth-size relation constitutes an age sequence. Lastly, wecompare our models to observations of giant molecular clouds and associated young star clusters inthe LMC and find good agreement between our model clouds and the observed relationship betweenH ii regions, young star clusters, and giant molecular clouds. Subject headings: stars: formation — ISM: clouds — ISM: evolution — turbulence INTRODUCTION
Giant molecular clouds (GMCs) are the primary reser-voir of molecular gas in the galaxy (Williams & McKee1997; Rosolowsky 2005; Stark & Lee 2006). Sincethe surface density of star formation shows a strongcorrelation with the surface density of molecular gas(Wong & Blitz 2002; Kennicutt et al. 2007; Bigiel et al.2008; Schruba et al. 2011), GMCs must also be the pri-mary site of star formation in the Milky Way. However,recent high-resolution observations have shown that theKennicutt-Schmidt law breaks down when the resolutionof an observation is finer than the typical length scales ofGMCs (Onodera et al. 2010; Schruba et al. 2010). Thus,in order to develop a detailed theoretical understandingof the relationship between star formation and molecu-lar gas, it is necessary to first understand the formation,evolution, and destruction of giant molecular clouds.One stumbling block in this effort is the substan-tial disagreement in the literature regarding both theformation mechanism and typical lifetimes of GMCs(see e.g. Goldreich & Kwan 1974; Zuckerman & Evans1974; Blitz & Shu 1980; Ballesteros-Paredes et al. 1999;McKee & Ostriker 2007; Murray 2011, and referencestherein). Some authors suggest that GMCs form outof bound atomic gas as a result of gravitational instabil-ity (Kim et al. 2002, 2003; Kim & Ostriker 2006; Li et al.2006a; Tasker & Tan 2009), surviving as roughly virial-ized objects for many cloud dynamical times (Tan et al.2006; Tamburro et al. 2008). In support of this pic-ture is the observation that massive clouds are found Department of Astronomy, 201 Interdisciplinary SciencesBuilding, University of California, Santa Cruz, CA, 95064, USA;[email protected] Department of Astronomy and Astrophysics, University ofToronto, Toronto, ON M5S 3H8, Canada Physics Department and Astronomy Department, Universityof California at Berkeley, Berkeley, CA 94720 to be marginally bound, with typical virial param-eters of order unity (Heyer et al. 2001; Rosolowsky2007; Roman-Duval et al. 2010). Since supersonicisothermal turbulence is found to decay via radiativeshocks in one or two crossing times (Mac Low et al.1998; Stone et al. 1998; Mac Low & Klessen 2004;Elmegreen & Scalo 2004), this model must invoke amechanism to drive supersonic motions for the life-time of a cloud, which could be several crossing times.Possible turbulent driving mechanisms include proto-stellar outflows (Norman & Silk 1980; McKee 1989;Li & Nakamura 2006; Li et al. 2010), H ii regions(Matzner 2002; Krumholz & Matzner 2009, hereafterKM09) , supernovae (Mac Low & Klessen 2004), or, asinvestigated here, mass accretion (Klessen & Hennebelle2010; V´azquez-Semadeni et al. 2010).Accretion driven turbulence in molecular clouds hasreceived little attention in the literature. However, asKlessen & Hennebelle (2010) point out, the kinetic en-ergy of accreted material can power the turbulent mo-tions observed in molecular clouds with energy conver-sion efficiencies of only a few percent. While there hasbeen no systematic study of the kinetic energy budgetof a molecular cloud formed via gravitational instability,this problem has been examined in the context of the for-mation of protogalaxies at high redshift. In one example,Wise & Abel (2007) analyzed simulations of virializinghigh redshift minihalos, tracking the thermal, kinetic,and gravitational potential energy of gas in protogalac-tic dark matter halos. In their models, which includeda nonequilibrium cooling model, gas collapsed onto thehalo and cooled quickly, causing turbulent velocities tobecome supersonic. As pointed out by Wang & Abel(2008), this means that the virialization process is a lo-cal one: gravitational potential energy can be converteddirectly into supersonically turbulent motions character-ized by a volume-filling network of shocks. The turbu- Goldbaum et al.lence in turn provides much of the kinetic support for thenewly virialized gaseous component of the dark matterhalo. If a similar mechanism is at work as gas cools andcollapses onto GMCs then gravitational potential energyalone may be sufficient to power turbulence in GMC.The most detailed simulations of simultaneously ac-creting and star forming giant molecular clouds wererecently completed by V´azquez-Semadeni et al. (2010).These numerical models included a simplified subgridprescription for stellar feedback by the ionizing radia-tion of newborn star clusters and focused on the bal-ance between accretion and feedback in clouds formedvia thermal instability in colliding flows. Throughout thecourse of the simulation, dense molecular gas condensedout of a warm atomic envelope, allowing a study of theinterplay between accretion and feedback in the simu-lated clouds. The resulting clouds were able to attain astate of quasi-virial equilibrium, in which the supply ofgas from the ambient medium balanced the formation ofstars and ejection of gas from H ii regions. Due to theidealized nature of the subgrid star formation feedbackprescription, in which all H ii regions were powered by acluster with the same ionizing luminosity, star formationfeedback was unable to act on the cloud as a whole butcould reduce the global star formation rate by destroy-ing overdensities. Since the simulation did not includestar clusters with large ionizing luminosities, the cloudas a whole could not be destroyed and star formationwould have eventually consumed all of the gas had thesimulation not been cut off. Even though the simulationsemployed a highly idealized star formation prescription,the computations still required substantial resources tocomplete and only allowed insight into the evolution ofa single cloud. It seems that a computationally inexpen-sive model that includes a somewhat more sophisticatedtreatment of star formation feedback is called for.In this work, we model the global evolution of giantmolecular clouds from their birth as low-mass seed cloudsto their dispersal after a phase of massive star formation.This is done using an updated version of the semianalyt-ical model of Krumholz, Matzner & McKee (2006, here-after Paper I). Using a virial formalism, we compute theglobal dynamical evolution of a single cloud while simul-taneously tracking its energy budget. Model clouds formstars, launch H ii regions and undergo accretion fromtheir environments. With these models, we are able toinvestigate the role accretion plays in maintaining turbu-lence in molecular clouds and directly compare to obser-vations of GMCs in the Milky Way and nearby externalgalaxies. This work is complementary to the simulationsof V´azquez-Semadeni et al. (2010), since our simplifiedglobal models allow us to survey a large variety of GMCsat little computational cost while including a much moresophisticated star formation feedback prescription. Weare able to capture model clouds with masses compara-ble to the most massive clouds observed in the MilkyWay and nearby galaxies, allowing us to simulate thesites of the majority of star formation in these systems(Williams & McKee 1997; Fukui & Kawamura 2010).We proceed by describing the formulation and imple-mentation of our GMC model in §
2. Next, in §
3, wetest our implementation of accretion. Following this, in § § §
6, we discuss the limitations inherent inthe simplifying assumptions we make to derive the cloudevolution equations. GOVERNING EQUATIONS
The GMC evolution model described below allows usto solve for the time-evolution of the global properties ofmodel molecular clouds. In contrast with previous work,we follow the flow of gas as it condenses out of the dif-fuse gas in the envelope surrounding the GMC and fallsonto the cloud. Employing simplifying assumptions aswell as the results of simulations of compressible MHDturbulence, we derive a set of coupled ordinary differ-ential equations that govern the time evolution of thecloud’s mass, radius, and velocity dispersion. Combin-ing the governing evolution equations with a set of ini-tial conditions, model parameters, and a model for thetime dependence of the mass accretion rate based on thegravitational collapse of the GMC envelope, we can solvefor the time evolution of the cloud. Below, we give anoverview of the model, discuss the formulation of our nu-merical scheme, describe our parameter choices, and givea brief description of our treatment of star formation andour model for the GMCs gas supply.
Model Overview
The model we employ below is based on the globalGMC model of Paper I, itself a generalization ofthe the global model for low mass star formation ofMcKee (1989), the Eulerian virial theorem (EVT) ofMcKee & Zweibel (1992), and the model for star formingclumps of Matzner (2001). Employing a virial formalism,we account for the dynamics and energy budget of gascontained within an Eulerian volume, V vir . We separatethe gas within V vir into three species: virial material, agaseous reservoir, and a photoionized wind. A schematicrepresentation of the components of our model is pre-sented in Figure 1.By design, each of the three components has a straight-forward physical interpretation. The first component,which we label virial material, consists of two physi-cally distinct subcomponents: a molecular cloud anda warm atomic envelope that encloses the cloud. Thecloud is assumed to be cold ( ∼
10 K), molecular, andcontained within a spherical volume of radius R cl . Theambient medium is composed of warm ( ∼ K) anddiffuse atomic gas that encloses the cloud and extendsbeyond the virial volume. The second component is agaseous reservoir, which we assume is composed of cold( ∼ K) neutral material that flows onto the cloud atfree-fall from beyond the virial volume. The last compo-nent is an ionized wind made up of hot ( ∼ K) ionizedgas ejected from the ionization fronts of blister-type H ii regions. All three components are allowed to mutuallyinterpenetrate. We restrict interaction between the com-ponents to the transfer of mass between the accretionflow and cloud as well as between the cloud and wind.Since the envelope and cloud are not allowed to inter-penetrate, we formally group the envelope and cloud to-gether; this somewhat artificial choice significantly sim-plifies the virial analysis.lobal Evolution of GMCs: Accretion 3 Figure 1.
A schematic overview of the GMC model. A molec-ular cloud ( black ) is embedded in a warm atomic envelope ( darkblue ). Cool atomic gas ( light blue ) flows onto the cloud, where itcondenses, recombines into molecules, and mixes with the cloud.Newborn OB associations ( blue stars ) drive H ii regions ( orange )and eject winds ionized winds back into the ambient medium. We make use of two simplifying assumptions regardingthe distribution and flow of the virial material. First,we assume that the virial material follows a sphericallysymmetric, smoothly varying density profile. Second, weassume that the the cloud is homologous: the cloud ex-pands, contracts, accretes, and sheds mass in such a wayas to always maintain the same smooth density profile.We assume a density profile of the form ρ ( r ) = ρ (cid:18) rR cl (cid:19) − k ρ for r ≤ R cl (1)where ρ is the density at the edge of the cloud and k ρ is assumed to be unity. This choice is consistent withthe Larson scaling relations observed in galactic (Larson1981; Solomon et al. 1987; Heyer et al. 2001, 2009) andextragalactic (Mizuno et al. 2001; Engargiola et al. 2003;Rosolowsky 2007; Bolatto et al. 2008; Hughes et al.2010) GMCs. Our assumed cloud density profile ef-fectively averages over clumpy and filamentary inter-nal structure and oblate shapes of observed clouds. At r = R cl , the density of the ambient medium is assumedto smoothly transition from ρ to ρ amb in a thin bound-ary layer. We assume that the density of the ambientmedium is negligible compared to the density of gasin the cloud, ρ amb ≪ ρ . The density of the ambientmedium remains ρ amb out to the virial radius.Beyond assuming a density profile, we must also spec-ify the velocity structure of all three components of themodel. We follow Paper I in assuming that the velocity ofthe virial material can be decomposed into a systematicand turbulent component, v = ˙ R cl R cl r + v turb . (2)We assume that v turb is randomly oriented with respectto position so that turbulent motions carry no net flux of matter. We make a similar assumption regarding thevelocity structure of the reservoir, v res = v res , sys ˆ r + v res , turb . (3)The systematic component of v res is due to the gravita-tional attraction of material within the virial volume v , sys Z r ∞ g · d r , (4)while the random component is such that( M − R V cl ρ v , turb dV ) / = √ σ res . Here g is thegravitational acceleration and σ res is the velocity disper-sion of gas in the reservoir feeding the accretion flow.Since the amount of material in the accretion flow isdetermined by how fast it can fall into the cloud, wemust simultaneously determine both the density profileand radial velocity of material in the accretion flow (seeAppendix B). Finally, for the wind material, we followPaper I in assuming v w = v + v ′ ej (5)where v ′ ej = 2 c ii ˆ r and c ii is the ionized gas sound speed.We follow McKee & Williams (1997) in choosing c ii =9 . − . Momentum Equation
In Appendix A, we derive the EVT for a simultane-ously evaporating and accreting cloud,12 ¨ I cl = 2( T − T ) + B + W − ddt Z S vir ( ρ v r ) · d S + a I ˙ M cl R cl ˙ R cl + 12 a I ¨ M cl R + a I ˙ M ej R cl ˙ R cl + 3 − k ρ − k ρ R cl ( ˙ M ej v ′ ej − ξ ˙ M acc v esc ) (6)Here, a I is a constant of order unity that depends onthe distribution of material in the cloud, I cl is the cloudmoment of inertia, T is the combined turbulent andthermal kinetic energy of the cloud, T is the energyassociated with interstellar pressure at the cloud sur-face, B is the net magnetic energy due to the pres-ence of the cloud, W is the gravitational term (equalto the gravitational binding energy in the absence ofan external potential [McKee & Zweibel 1992]), the sur-face integral is proportional to the rate of change ofthe moment of inertia inside the bounding viral surface, M cl is the cloud mass, R cl is the cloud radius, ˙ M ej isthe mass ejection rate, ˙ M acc is the mass accretion rate, v esc = { G [ M cl + M res ( R cl )] /R cl } / is the escape ve-locity at the edge of the cloud, and ξ is a dimensionlessfactor we compute via equation (A8) that depends on thedepth of the cloud potential well. The quantity ξv esc isthe accretion rate weighted average infall velocity. Pre-cise definitions for T , T , B , and W are given in PaperI. The EVT of a cloud without accretion or mass losswould only contain the terms up to the surface integral.The next three terms account for changes in the cloudmoment of inertia due to changes in the mass of thecloud, while the last term accounts for the rate at which Goldbaum et al.the recoil of inflowing and outflowing mass injects mo-mentum into the cloud. Inflows and outflows are treatedseparately because material is ejected at a constant ve-locity, but is accreted at a velocity that is a function ofthe depth of the potential well of the cloud. The dimen-sionless factor ξ appears due to this difference.The mass of the cloud can only change via mass accre-tion or ejection, ˙ M cl = ˙ M ej + ˙ M acc . (7)We assume that ejection of material can only decreasethe mass of the cloud and accretion can only increasethe mass of the cloud. Since stars may not follow thehomologous density profile we assume, we neglect thechange in the cloud mass due to star formation. Weexpect the error incurred from this assumption will besmall, since stars make up a small fraction of the mass ofobserved clouds (Evans et al. 2009; Lada et al. 2010) andour assumed star formation law converts a small fractionof the cloud’s mass into stars per free-fall time.We follow Paper I in using the assumption of homol-ogy and the results of simulations of MHD turbulence toevaluate each term in the EVT in terms of constants andthe dynamical variables R cl , M cl , and σ cl . In the end,we obtain a second order nonlinear ordinary differentialequation in R cl , a I ¨ R cl = 3 c R cl + 3 . σ R cl − a ′ (1 − η B ) GM cl R − πP amb R M cl − a I ˙ M acc M cl ˙ R cl + (cid:18) − k ρ − k ρ (cid:19) ˙ M ej M cl v ′ ej − ξ ˙ M acc M cl v esc ! . (8)Here, a ′ = 15 − k ρ − k ρ (cid:20) − k ρ ) Z x − k ρ y ( x ) dx (cid:21) , (9) y ( x ) is the ratio of the mass of reservoir material to themass of cloud material contained within a normalizedradius x = r/R cl (see Appendix B), and η B is the ratioof the magnetic critical mass to the cloud mass.This equation governs the balance of forces acting onthe cloud as a whole. Each term corresponds to a sin-gle physical mechanism that can alter the radial forcebalance. The first two terms are due to thermal and tur-bulent pressure support, respectively. The third is due toa combination of gravitational compression and magneticsupport. The fourth is due to the confining interstellarpressure. The fifth comes from the exchange of momen-tum between the expanding cloud and infalling accretionflow. Finally, the last term is due to a combination ofthe recoil from ejected material and the ram pressure ofaccreting material. Although the two parts of the recoilterm have opposite signs, ˙ M ej and ˙ M acc have oppositesigns as well: ˙ M ej < M acc >
0. This implies thatboth the recoil due to launching wind material and theram pressure of accreting reservoir material tend to con-fine the cloud.Letting M cl , , R cl , , and σ cl , be the cloud mass, radius,and velocity dispersion at t = 0 and defining the initial cloud crossing time, t cr , = R cl , /σ cl , , we can define thedimensionless variables M = M cl /M cl , , R = R cl /R cl , , σ = σ cl /σ cl , , and τ = t/t cr , . Letting primes denote dif-ferentiation with respect to τ , we can write equation (8)in dimensionless form R ′′ = 3 . σ + 3 M − a I R − η G a ′ MR − η P R M − M ′ acc R ′ M + η E M ′ ej M − η A ξM ′ acc ( f M R ) / (10)where M = σ cl , /c cl (11)is the initial turbulent Mach number and we define thedimensionless constants η G = 3(1 − η B ) a I α vir , (12) η P = 4 πR , P amb a I M cl , σ , (13) η E = (cid:18) − k ρ − k ρ (cid:19) v ′ ej σ cl , (14) η A = (cid:18) − k ρ − k ρ (cid:19) (cid:18) α vir , (cid:19) / (15)where α vir , = 5 σ , R cl , GM cl , (16)is the initial nonthermal virial parameter(Bertoldi & McKee 1992). These constants are setby the ratios of various forces acting on the initialstate of the cloud. η G is proportional to the ratio ofthe initial magnetic forces to the initial gravitationalforce, and η P , η E , and η A are the ratios of the ambientpressure force, the mass ejection recoil force, and theinitial accretion ram pressure force to the initial internalturbulent forces, respectively.Comparing equation (10) with the corresponding equa-tion given in Paper I, we see that two new terms pro-portional to M ′ acc have appeared. In practice, we findthat, of the two terms, the one proportional to η A domi-nates, implying that the primary direct impact of accre-tion on the radial force balance of the cloud is to providea confining ram pressure. We will see in the next sec-tion that accretion also increases the turbulent velocitydispersion, implying that the kinetic pressure term alsoincreases when accretion is included. The cloud radius isdetermined by a balance between kinetic pressure and acombination of gravity, accretion ram pressure, and windrecoil pressure. Thermal pressure support is negligible.We also note that although the ambient pressure termis of the same form as in Paper I, we assume an ob-servationally motivated value for the ambient pressure, P amb /k B = 3 × K cm − (McKee 1999). This includesthermal and turbulent pressure but neglects magneticand cosmic ray pressure, since magnetic fields and cos-mic rays permeate both the cloud and the ambient ISM.We have also adjusted the ambient pressure upwards bya factor of two because GMCs form in overdense regionsof the ISM where the hydrostatic pressure is higher thanlobal Evolution of GMCs: Accretion 5average. In Paper I, P amb was chosen to be artificiallyhigh to ensure that the cloud would start its evolutionin hydrostatic equilibrium. This choice was made to ac-count for the weight of the gaseous reservoir that wasnot explicitly included. In practice, by choosing a lowervalue for P amb , we find that the ambient pressure termis subdominant for most of the evolution of the cloud.This is expected, since we now correctly account for thepressure of the reservoir through the term proportionalto η A . Once accretion halts, the cloud is left out of pres-sure equilibrium and must expand to match the ambientpressure. This effect is seen most clearly in the secondcolumn of Figure 3. Energy Equation
In Appendix C, we derive the time evolution equationfor the total energy of the cloud, d E cl dt = ˙ M cl M cl [ E cl + (1 − η B ) W ] − πP amb R ˙ R cl + GM cl ˙ M cl R cl χ − M cl ˙ R cl ˙ M cl R cl ! + (cid:18) − k ρ − k ρ (cid:19) ˙ R cl ( ˙ M ej v ′ ej − ξ ˙ M acc v esc )+ ϕ ( 32 ˙ M acc σ + 32 ˙ M acc σ + γ ˙ M acc v ) − a I ˙ M acc ˙ R − M acc σ cl + G cl − L cl , (17)where G cl and L cl are the rates of energy gain and lossdue to H ii regions and turbulent dissipation, respec-tively, E cl is the total energy due to the presence of thecloud (see equation [C5]), σ res is the velocity dispersionof material that is being accreted, χ is given by Equa-tion C8, γ is given by equation (C15), and ϕ is a freeparameter that sets the amount of energy available todrive accretion driven turbulence. The parameter ϕ isthe only adjustable constant in our model that is notconstrained by the results of simulations or observationsand must be tuned to reproduce the observed propertiesof clouds. The evolution of the cloud is very sensitive to ϕ and we justify our fiducial choice, ϕ = 0 .
75, in § M cl /M cl ) E cl is the rate of change ofthe cloud’s energy as mass is advectively added to or car-ried away from the cloud. Similarly, ( ˙ M cl /M cl )(1 − η ) W is the rate of change of the gravitational and magneticenergy due to changes in the mass of the cloud. The nextterm is the rate at which external pressure does compres-sional work on the cloud. This is followed by a term thataccounts for the gravitational work done on the cloud bythe reservoir as the cloud expands and contracts. Thefollowing term represents the rate at which mass inflows,outflows, and external thermal and turbulent pressure,respectively, do compressional work on the cloud. Thenext term, which is proportional to ϕ , represents the rateof kinetic energy injection via stirring of turbulence byaccreted material. This is followed by two terms that areproportional to ˙ M acc , which account for the fact that inthe frame comoving with the motions of material in thecloud, accreted material is moving at the transformed velocity, v res − v , different than the velocity of the reser-voir material in the rest frame, v res . Lastly, G cl and L cl are the rate of energy injection by star formation and therate at which energy is radiated away, respectively.Noting that turbulent motions carry no net radial fluxof matter and recalling that we had set B turb = 0 . T turb ,we may evaluate equation (C5) and obtain for the totalcloud energy, E cl = 12 a I M cl ˙ R + 2 . M cl σ + 32 M cl c − (cid:20) a ′ (1 − η B ) + χ (cid:21) GM R cl . (18)Taking the time derivative of this expression, substitut-ing into equation (17), and nondimensionalizing as in § σ ,4 . a I σ ′ = − R ′ R ′′ σ − η G M R ′ R σ − η P R R ′ M σ + η E M ′ ej R ′ M σ − η A M ′ acc R ′ ( M R ) / σ − M ′ acc R ′ M σ − (3 − . ϕ ) M ′ acc σa I M + η D ϕςM ′ acc M σ + η I ϕγM ′ acc f Rσ + G − L a I M σ (19)where ς = σ res /σ res , , G − L = R cl , ( G cl − L cl ) M cl , σ , , (20)and we define the constants, η D = 3 σ , a I σ , , (21) η I = 10 a I α vir , . (22)Here, η D is propotional to the ratio of the initial turbu-lent kinetic energy in the reservoir and the initial turbu-lent kinetic energy in the cloud and η I is proportional tothe ratio of the initial kinetic energy due to the infall ofthe reservoir to the initial turbulent kinetic energy of thecloud.Since motions in GMCs are highly supersonic, theinternal structure of a typical cloud is characterizedby strong shocks. Because clouds have short coolingtimescales, the shocks present throughout GMCs must beradiative. The braking of turbulent motions via radiativeshocks has been extensively studied in numerical simula-tions (see e.g. Mac Low et al. 1998; Stone et al. 1998) inwhich the turbulent dissipation timescale is found to be t dis = E turb / ˙ E = kλ in /σ cl where k is a constant of orderunity and λ in is the characteristic length scale of tur-bulent energy injection. The simulations of Stone et al.(1998) give k = 0 .
48 and E turb = 2 . M cl σ . Motivatedby this result and using a scaling argument given byMatzner (2002) and McKee (1989), we assume that thedimensionless rate of energy loss is given by L = η v φ in M σ R . (23) Goldbaum et al.Here η v is a constant of order unity that depends on thenature of MHD turbulence in the cloud and we assume φ in = λ in / R cl . The factor of 4 in our expression for φ in comes from the fact that the largest wavelength modesupported by the cloud is λ max = 4 R cl , corresponding tonet expansion or compression of the cloud. We make useof the simulations of Stone et al. (1998) to calibrate thisexpression. For the run most resembling real molecularclouds, we find η v = 1 .
2. We follow Paper I in adopting φ in = 1 . λ in cannot be directlyobserved, with the velocity structure of simulated clouds,where λ in is known a priori , and found λ in & R cl .Comparing our velocity dispersion evolution equation(equation [19]) to the corresponding equation given inPaper I, we see there are four new terms proportional to M ′ acc . In practice, we find that the primary effect of ac-cretion on the energy balance of the cloud is to increasethe turbulent velocity dispersion via the terms propor-tional to ϕ . We will show in § Star Formation and H ii Regions
Star formation is able to influence the evolution of thecloud by ejecting mass and by injecting turbulent kineticenergy as expanding H ii regions merge with and driveturbulent motions in the cloud. The first mechanismis accounted for in our models by including an ionizedwind that decreases the mass of the cloud and confinesthe cloud by supplying recoil pressure. The second mech-anism is accounted for by the G cl term in equation (17)that represents energy injection by H ii regions.Since we only know the global properties of thecloud, we calculate the rate of star formation by mak-ing use of a power-law fit to the star formation lawof Krumholz & McKee (2005). Stars form at a lowefficiency per free-fall time, consistent with observa-tions of star formation in nearby molecular clouds(Krumholz & Tan 2007). Individual star formationevents occur once a sufficient amount of mass has ac-cumulated to form a star cluster. Star cluster massesare found by drawing from a cluster mass function ap-propriate for a single cloud (see equation [44] of PaperI). We then populate the cluster with individual stars bypicking masses from a Kroupa (2002) IMF. If the totalionizing luminosity of the newborn star cluster is suffi-cient to drive the expansion of an H ii region, we beginto track the resulting expansion.Once a massive star cluster forms, it photoionizes gasin its surroundings and drives the expansion of an H ii region. Paper I tracked the expansion of individual H ii regions by assuming the analytic self-similar solution forH ii region expansion worked out by Matzner (2002).This solution uses the fact that once an H ii region hasexpanded beyond the Str¨omgren radius, most of the massin the H ii region volume is in a thin shell of atomic gasat a radius r sh from the center of the H ii region. Theionized gas in the interior of the shell exerts a pressureon the surface of the shall, causing the shell to accelerateoutwards. The shell evolution equation derived from thisanalysis admits a self-similar solution for the expansion of the H ii region. This self-similar result does a good jobof predicting the expansion if there is no characteristicscale in the problem.However, the introduction of radiation pressure leadsto a characteristic radius, r ch , and time, t ch , at which thegas pressure and radiation pressure at the inner surfaceof the shell are equal. Radiation pressure is the dominantforce driving the expansion of the ionized bubble when r sh < r ch and gas pressure dominates when r sh > r ch .KM09 modified the theory of Matzner (2002) to accountfor the effect of radiation pressure in the initial stages ofthe expansion. They derived an explicit functional formfor r ch and t ch in terms of the bolometric and ionizingluminosity of the central star cluster, properties of themolecular cloud, and fundamental constants (see equa-tions [4] and [9] in KM09). The numerical value of r ch and t ch depends on the the bolometric luminosity of thecentral star cluster and the ionizing photon flux of thecentral star cluster. The value we choose for f trap , afactor that accounts for the trapping and reradiation ofphotons as well as the trapping of main sequence windswithin the neutral shell, and φ , a factor that accountsfor the absorption of radiation by dust, are the fiducialvalues quoted by KM09.Defining the dimensionless variables x sh = r sh /r ch and τ sh = t/t ch , the equation of motion for the shell reducesto (KM09) ddτ sh (cid:18) x dx sh dτ sh (cid:19) = 1 + x / . (24)This assumes that gas in the neighborhood of the H ii region follows a density profile proportional to r − , effec-tively placing the H ii region in the center of the cloud.This accounts for the fact that H ii regions form in over-dense regions of the cloud. In practice, we solve equa-tion (24) numerically to obtain x sh and dx sh /dτ sh andthus r sh ( t ) and ˙ r sh ( t ).In a gas pressure driven H ii region, the force exertedon the expanding bubble is twice the recoil force pho-toionized material imparts on the cloud as it is ejected(Matzner 2002). The photoevaporation rate can then bestraightforwardly calculated via ˙ M ej = − ˙ p sh / c ii . Here p sh is the momentum of the shell and ˙ p sh is the forceacting on the shell. In a radiation pressure driven H ii region, the total force is given by the sum of the gas pres-sure and radiation pressure forces, ˙ p sh = ˙ p gas + ˙ p rad . Atearly times, when r sh ≪ r ch , the radiation force domi-nates, so ˙ p rad ≫ ˙ p gas . Thus, If we calculate the massejection rate from the the total force acting on the shell,we will overestimate the mass ejection rate in a radiationpressure dominated H ii region. To correct for this effect,we modify the analysis of Paper I by only including thegas pressure force when we calculate the mass ejectionrate, ˙ M ej = − ˙ p gas c ii (25)lobal Evolution of GMCs: Accretion 7where, ˙ p gas = (cid:18) Sφπα B (cid:19) / . k B T ii r / (26)= 3 . × L x / dynes . (27)Since we do not include the dynamical ejection of ma-terial as H ii regions break out of the cloud surface,this is formally a lower limit on the true mass ejectionrate. Since clouds are clumpy and somewhat porous(Lopez et al. 2011), we expect to make little error byneglecting dynamical ejection.Once the stars providing 50% of the total ionizing lu-minosity of the central star cluster have left the mainsequence, the H ii region enters an undriven momentum-conserving snowplow phase. When the expansion veloc-ity of the H ii region is comparable to the turbulent ve-locity dispersion, we assume that the H ii region breaksup and contributes turbulent kinetic energy to the cloud.H ii regions can merge with the cloud during either thedriven or undriven phases. If the radius of the H ii regionis greater than the radius of the cloud and the expansionvelocity of the shell is greater than the cloud escape ve-locity, we say the H ii region disrupts the cloud and endthe global evolution.If an H ii region merges with the cloud at time t = t m when its radius is r sh = r m , the rate of energy injectionfrom a single H ii region is given by G cl = 1 . η E T ( t m ) (cid:18) r m R cl (cid:19) / δ ( t − t m ) . (28)Here T ( t m ) = p sh σ cl / η E paramaterizes the efficiencyof energy injection, δ ( t ) is the Dirac delta function,and the factor of 1.6 arises because magnetic turbu-lence is slightly sub-equipartition compared to kineticturbulence. The factor ( r m /R cl ) / accounts for themore rapid decay of turbulence when the driving scaleis smaller (see Paper I). Mass Accretion
Consistent with the analysis in Appendix B, we treatthe reservoir as a gravitationally unstable spherical cloudundergoing collapse. The cloud is primarily composed ofatomic gas in both warm and cold phases. We expectthat as the reservoir collapses, material that is accretedonto the cloud will cool and become molecular. We ap-proximate that the reservoir has approximately constantsurface density, Σ res . To find an upper limit on the massaccretion rate, we can assume that the reservoir is under-going collapse in the limit of zero pressure. It is straight-forward to show that the resulting accretion rate ontothe central condensation is given by˙ M acc , ff = 256 π G Σ t (29)Since the estimate for the accretion rate given in equa-tion 29 does not take into account pressure support,it is likely an overestimate of the true accretion rate.Indeed, McKee & Tan (2003) considered the inside outcollapse of equilibrium polytropic molecular cloud coresand, in the case of k ρ = 1, found the same scaling with time and surface density but a substantially lower coeffi-cient. Tan & McKee (2004) argued that subsonic inflowwas a more realistic initial condition condition than thestatic equilibrium solution used by McKee & Tan (2003).Hunter (1977) found that the set of solutions for the col-lapse of an isothermal sphere starting at an infinite timein the past. The solution with an infall velocity of about c cl / . k ρ = 1corresponds to γ = 0 (McKee & Tan 2003) rather thanHunter’s γ = 1. Nonetheless, we assume that the accre-tion rate for our problem is also 2.6 times greater thanthat for a static initial condition and find˙ M acc , TM04 = 10 . G Σ t . (30)Clearly, this result is uncertain, and magnetic fields in-troduce further uncertainty. Fortunately, we find thatvarying the numerical coefficient in equation 30 does notaffect the qualitative nature of the results discussed be-low.To model the effect of a finite gas supply, once thetotal mass of gas that has fallen onto the cloud ex-ceeds the total mass of the reservoir, M res , we set˙ M acc = 0. When comparing with galactic populationsof GMCs, we set M res = 6 × M ⊙ , the observedupper mass cutoff for GMCs (Williams & McKee 1997;Fukui & Kawamura 2010). This may underestimate thetrue upper mass since fragmentation may lead to a rangeof reservoir masses. Since 6 × M ⊙ GMCs are rela-tively rare, we set M res = 2 × M ⊙ for the runs pre-sented in Figure 3 and Table 2 as ∼ M ⊙ is a moretypical GMC mass.If we also assume that the atomic reservoir is virializedsuch that the virial parameter of the reservoir, α res , isconstant with radius, we find σ res = 0 . α / G Σ res t (31)Fiducially, we take α res = 2 .
0, corresponding to amarginally gravitationally bound reservoir. This param-eterization assumes that the reservoir satisfies an internallinewidth-size relation of the form σ res ( r ) ∝ r / . Theincrease of the reservoir velocity dispersion with time re-flects that material originates at increasingly larger radii.The precise normalization of the mass accretion law isa major source of uncertainty in our modeling. Since thelength scales over which material is swept up into thecloud through the reservoir approach galactic dynamicalscales (Dobbs et al. 2011; Tasker 2011), a more completetreatment would require tracking the gas dynamics fromthe scale of an entire galaxy down to the the scale of thereservoir. This would also allow us to self-consistentlymodel the end of accretion onto the cloud instead of as-suming that mass accretion cuts off abruptly. In a forth-coming paper, we plan to add our molecular cloud mod-els to a simulation of gas dynamics in a galactic disk tomodel the gas reservoir for a population of GMCs.Although the normalization of the accretion law issomewhat uncertain, we can make use of observations Goldbaum et al.of the gas content of nearby galaxies to estimate Σ res ,the surface density of gas in the reservoir feeding thecloud. The average H i surface density in the inner diskof the Milky Way is observed to be approximately con-stant, ∼ ⊙ pc − (Kalberla & Dedes 2008). Beyonda galactocentric radius of ∼
10 – 15 kpc, the H i surfacedensity exponentially decreases with radius, however, fewGMCs are observed in the outer Milky Way (Heyer et al.2001) or beyond the optical radius of nearby galaxies(Engargiola et al. 2003; Bigiel et al. 2008). Similar satu-rated mean H i surface densities are observed in nearbygalaxies, except in the central regions of some galaxieswhere the ISM becomes fully molecular and the H i sur-face density goes to zero (Leroy et al. 2008). For thisreason, we adopt Σ res = 8 M ⊙ pc − as a fiducial atomicreservoir surface density typical of the bulk ISM of localstar forming galaxies.Although the mean atomic surface density in nearbygalaxies is as low as 8 M ⊙ pc − , the atomic ISM is ob-served to be clumpy, with overdense regions reaching sig-nificantly higher surface densities. These regions may beassociated with spiral arms, as in M 33 (Thilker et al.2002), or driven by gravitational instability, as in theLMC (Yang et al. 2007). For this reason, we also explorethe behavior of molecular clouds accreting from highersurface density gas, Σ res = 16 M ⊙ pc − . Since massivemolecular clouds are universally observed to be associ-ated with high surface density gas (Wong et al. 2009;Imara et al. 2011), we expect there to be marked differ-ences between molecular clouds that accrete from highsurface density gas and clouds that accrete from lowsurface density gas. Although the gas will be primar-ily atomic at a gas surface density of 16 M ⊙ pc − , weexpect that there should be some diffuse “dark” molecu-lar gas (Krumholz et al. 2008; Wolfire et al. 2010). Thusthe reservoir is not necessarily completely atomic, butinstead primarily composed of atomic gas. Numerical Scheme
Equations (7), (10), and (19) constitute a system ofcoupled, stochastic, nonlinear ordinary differential equa-tions in M , σ , and R . We solve these equations by usinga straightforward Euler integration with an adaptive stepsize. The precise order in which we update cloud proper-ties is as follows. After calculating the instantaneous starformation rate, we calculate M ′ by summing the compo-nents due to ionized winds and mass accretion. Next,we calculate the rate of turbulent dissipation using equa-tion (23). We then calculate ζ , the ratio of the massdoubling time to the free-fall time, using equation (B5).We use ζ to calculate a ′ , f , and ξ by interpolating onprecomputed tables. Since σ ′ depends on R ′′ , we firstevaluate R ′′ using equation (10) and then compute σ ′ using equation (19). Next, we check if R , M , or σ willchange by more that 0 . R , M , and σ are smaller than 0 . R , M , and σ at the new time step, update the state ofany H ii regions created in previous time steps, and thencreate new H ii regions using the procedure described in § • The time step is less than 10 − of the current evo-lution time (i.e. ∆ τ /τ < − ). • The mean visual extinction falls below A V, min =1 .
4, corresponding to the CO dissociation thresholdfound by van Dishoeck & Black (1988). • An H ii region envelops and unbinds the cloud, i.e.if r sh > R cl and ˙ r sh > v esc .We use the phrases collapse, dissociation, and disruption,respectively to describe these scenarios. The dissociationthreshold depends on the ambient radiation field. How-ever, Wolfire et al. (2010) found that the CO dissociationthreshold varies by only a factor of two when the intensityof the ambient radiation field varies by an order of mag-nitude, so we neglect variations in the radiation field andfor consistency with Paper I, adopt A V, min = 1 .
4. Thesurface density corresponding to the dissociation thresh-old depends on the assumed dust to gas ratio and thuson the metallicity. For solar metallicity, a surface densityof 1 g cm − corresponds to A V = 214 .
3. Since we use adissociation threshold based on CO rather than H todefine the end of the cloud’s life, we may miss the fur-ther evolution of a diffuse molecular cloud where most ofthe carbon is neutral or singly ionized but the hydrogenis still molecular.These halting conditions probably oversimplify thetrue end of a cloud’s evolution due to our assumptionof spherical symmetry and homology. For the case ofcollapse, it is more likely that the cloud would undergorunaway fragmentation rather than monolithic collapse.In the case of dissociation, even if the mean surface den-sity drops below the point where CO can no longer re-main molecular, that does not preclude the possibilitythat overdense clumps might retain significant amountsof CO. Finally, for the case of disruption, even if an H ii region delivered a large enough impulse to unbind thecloud, it may simply be displaced as a whole, or be dis-rupted into multiple pieces which would then evolve in-dependently. It is also likely that if the cloud is disruptedwhile still actively accreting, the cloud would inevitablyrecollapse since it is unlikely that an H ii region wouldhave enough kinetic energy to unbind the reservoir.Since we must necessarily use simple criteria to haltthe cloud evolution, our estimates of cloud lifetimes pre-sented below are strictly lower limits to the true lifetimeof a cloud. Both the disruption and dissociation criteriado not preclude the presence of overdense clumps thatmay survive the destruction events. However, our esti-mates of cloud lifetimes are appropriate for the lifetimeof a single monolithic cloud. Any overdense clumps thatdo survive would represent entirely different clouds thatwould evolve independently of each other. Input Parameters
To complete our model, we must choose a set of param-eters and initial conditions to fully determine our cloudevolution equations. In Table 1 we have listed the vari-ous fiducial parameters we have chosen for our model aswell as the references from which we derive our choices.Some of these parameters are motivated by observations,lobal Evolution of GMCs: Accretion 9
Table 1
Fiducial ParametersParameter Value Reference α vir , cl ,
60 M ⊙ pc − Heyer et al. (2009) c s a .
19 km s − — c ii b .
74 km s − McKee & Williams (1997) η B η E η v φ in ϕ A v , min . ⊙ pc − van Dishoeck & Black (1988) α res P amb /k B × K cm − McKee (1999) M cl , × M ⊙ — a Assumes T = 10 K and µ = 2 . b Assumes T = 7000 K and µ = 0 . others by the results of simulations, and one parameter( ϕ , see §
3) is left free.The cloud initial conditions can be computed givenan initial cloud mass along with our assumed value for α vir , and Σ cl , . Matzner & McKee (2000) and Fall et al.(2010) found that protostellar outflows are energeticallyimportant when M cl . . M ⊙ . Above this mass, theycontribute negligibly. Thus, if we choose a low initialmass, we may be underestimating the amount of turbu-lent energy injection by star formation feedback at earlytimes since we do not account for protostellar outflows.For this reason, we choose a relatively large initial mass, M cl , = 5 × M ⊙ . For reference, given our choiceof initial virial parameter and surface density, this cor-responds to, R cl , = 11 . σ cl , = 2 . − , and t cr , = 4 . MODELS WITH ACCRETION ONLY
Before beginning full simulations using our method, wemust first test the behavior of our model as we vary thefree parameter ϕ introduced in the derivation of the en-ergy evolution equation. For this purpose, we have runour model with star formation feedback disabled. Theonly physical mechanisms modeled in these tests are ac-cretion and the decay of turbulence. Since there is norandom drawing from the stellar or cluster IMF in theseruns, the results are fully deterministic. These simplifiedmodels allow us to understand how the results depend onour choice for the tunable parameter ϕ and provide phys-ical insight that will be useful in interpreting the resultsof the more complex and stochastic runs that includefeedback.The energetics and virial balance of our cloud modelsdepend critically on the parameter ϕ . Broadly speak-ing, ϕ controls the amount of turbulent kinetic energyinjected by the accretion flow. For the case ϕ = 0 the ac-cretion flow contributes the minimum possible amount ofturbulent kinetic energy. This means that accreted mate-rial cannot contribute significantly to turbulent pressure Figure 2.
Cloud surface densities ( bottom row ), virial parame-ters ( second row ), velocity dispersions ( third row ), and radii ( toprow ) for 400 different runs, each with a different choices for ϕ , asindicated in the color bar. Star formation was turned off for allruns. support, since the accreted material is maximally sub-virial. With this choice, as the cloud accretes mass, itsenergy budget must become increasingly dominated byself-gravity. Once this happens, internal pressure sup-port is negligible and the cloud must inevitably undergogravitational collapse.Alternatively, we could set ϕ >
0. If ϕ is small, theturbulent kinetic energy of a newly accreted parcel ofgas would still be small compared to the gravitationalpotential energy of the gas parcel. Thus, once the cloudis primarily composed accreted gas, the cloud will un-dergo gravitational collapse, although on a slightly longertimescale than in the ϕ = 0 case. At some larger value of ϕ , accretion contributes a net positive amount of energyto the cloud, balancing out the negative gravitational po-tential energy of the newly accreted material. The couldwill still collapse with this choice since turbulent motionsquickly decay away. As ϕ is increased further, we shouldeventually find that at some critical value, ϕ = ϕ crit , ac-cretion drives turbulent motions with sufficient vigor toavoid the gravitational collapse of the cloud entirely.The results of test runs with different choices of ϕ arepresented in Figure 2. The time evolution of the cloudsurface density, virial parameter, velocity dispersion, andcloud radius is plotted for a selection of clouds evolvedwith different choices for ϕ . Each line depicts the timeevolution of a cloud property and is color coded by thevalue of ϕ chosen. It is obvious that the value of ϕ canstrongly influence the resulting evolution.If ϕ = 0, the cloud experiences global collapse in afree-fall time. Initially, the cloud velocity dispersiondecreases, but inevitably the gravitational term in thevelocity dispersion evolution equation, proportional to − R ′ /R , becomes dominant, and the velocity dispersionbegins to diverge. The fact that σ cl diverges as R cl goesto zero is an artifact. In reality, the highest density re-gions would independently fragment and collapse and the0 Goldbaum et al.cloud would never undergo a monolithic collapse.As we increase ϕ , the cloud is able to support itselfagainst collapse for longer periods. Near ϕ = 0 .
5, accre-tion brings in net positive energy but turbulent dissipa-tion wins out, and the cloud still eventually collapses. Ata critical value, ϕ crit ≃ .
8, accretion driven turbulencealone is sufficient to hold up the cloud against collapse foras long as the reservoir continues to supply mass to thecloud. The mass, radius, and velocity dispersion of thecloud increase in such a way as to maintain a constantvirial parameter and surface density.Since we expect that gas motions driven by accret-ing dense clumps should be at least somewhat correlatedwith the motions of the infalling clumps, we do not ex-pect a physically realistic choice of ϕ to be very close tozero. On the other hand, a model in which a cloud is en-tirely supported by accretion driven turbulence seems topreclude the possibility that a significant fraction of thekinetic energy of infalling gas is radiated away in an ac-cretion shock. For this reason, we rule out as unphysicalruns with ϕ ≈ ϕ ≥ ϕ crit . The precise value of ϕ wewill use in our models that include star formation belowdepends on uncertain details of the accretion and mix-ing of infalling gas. In practice, we find that even withthe energy provided by star formation feedback, cloudsgenerally undergo free-fall collapse or reach unreasonablyhigh mean surface densities once they are primarily com-posed of accreted material if we choose ϕ . .
7. Sinceclouds are generally not observed to be in global free-fall collapse, we instead pick a value somewhat higherthat this, ϕ = 0 .
75 for our fiducial models. This splitsthe difference between accretion contributing a negligi-ble amount of energy to the cloud when ϕ = 0 . ϕ = 1. We will see below that our fidu-cial choice broadly reproduces the observed properties ofmolecular clouds in the Milky Way and nearby galaxies. MODELS WITH ACCRETION AND STAR FORMATION
Feedback by the action of ionizing radiation emittedby newborn stellar associations alters the evolution of aGMC after the birth of the first massive star cluster. Thesource of energy provided by massive star formation canbe a significant component of the energy budget of theentire cloud. For the remainder of this paper, we considermodels with the star formation prescription described in § Overview of Results
We have run two sets simulations with parameters cho-sen to model conditions in interarm (Σ res = 8 M ⊙ pc − )and spiral arm (Σ res = 16 M ⊙ pc − ) regions. Besidesthe two different choices for the ambient surface density,all other parameters and initial conditions are identical.The time evolution of a subsample of runs are plottedin Figure 3 and average properties of the full sample arepresented in Table 2.The most striking result of our comparison is that thefinal mass of our model molecular clouds depends on theassumed mass accretion history. Clouds evolved with alow accretion rate, corresponding to conditions in inter-arm regions, grow larger than 10 M ⊙ less than 30% ofthe time and very rarely reach masses comparable to the most massive GMCs in the local group. The vast major-ity of clouds are instead disrupted by an energetic H ii region within a few crossing times. The clouds attaina quasi-equilibrium configuration in which mass accre-tion is roughly balanced by mass ejection. Clouds avoidglobal collapse by extracting energy from the expansionof H ii regions.The evolution of the clouds is characterized by discreteenergy injection events due to the formation of a sin-gle massive star cluster. Once a cluster forms, it ejectsa wind and launches an H ii region. The recoil forceof launching the wind leads to an overall confining rampressure, causing the radius to decrease and the surfacedensity to increase. Once the star cluster burns out, theH ii region expansion decelerates and then stalls. Whenthe expansion velocity of the H ii region is comparableto the cloud velocity dispersion, the kinetic energy of theexpanding H ii region is converted into turbulent kineticenergy, causing a spike in the turbulent velocity disper-sion. The turbulent kinetic energy exponentially decaysaway over a crossing time, but the temporarily elevatedvelocity dispersion increases the turbulent kinetic pres-sure, causing the cloud to expand. This leads to oscilla-tions in the cloud radius and mean surface density. Onthe whole, clouds that are not quickly disrupted by H ii regions, are able to survive as quasi-virialized objects forseveral crossing times before they are either disrupted ordissociated.Clouds evolved with a higher ambient surface density,typical of spiral arm regions in the Milky Way, exhibitsignificantly different behavior. Since these clouds ac-crete mass much faster than in the low surface densityruns, they are not able to attain steady state between ac-cretion and ejection of mass. While some clouds are stilldestroyed by energetic H ii regions early in their evolu-tion, over 90% of these clouds were able to accrete theirentire reservoir after 25 Myr. At this point, the cloudsare generally quite massive, ∼ . × M ⊙ . Once ac-cretion is shut off, the clouds are no longer confined byaccretion ram pressure and lose a portion of the powerthat had been driving turbulence. For this reason, thevelocity dispersion decreases in response to the loss of ac-cretion driven turbulence, and the cloud radius expandsin response to the loss of the confining pressure providedby accretion. Before the cloud can dissociate, it attainspressure balance with the ambient ISM at a lower ve-locity dispersion and larger radius. For the next 20 to30 Myr, the clouds evolve in much the same way as themassive cloud models considered in Paper I. The cloudscan be supported against self-gravity for many dynam-ical times by forming stars and launching H ii regions.Particularly energetic H ii regions can disrupt the cloudsand excursions to low surface density can dissociate theclouds. The lifetime of these clouds is thus set by theamount of time they can accrete. This may imply thatspiral arm passage times set GMC lifetimes, althoughfurther work is needed to clarify this tentative conclu-sion.The star forming properties of the two sets of modelsare also somewhat different. In the interarm case, thestar formation efficiency, ǫ = M ∗ , tot R t life ˙ M acc dt , (32)lobal Evolution of GMCs: Accretion 11 Table 2
Average Properties of Model CloudsΣ res h t life i h M max i (cid:10) M ejected (cid:11) h ǫ i h ǫ ff i N dissoc N disrupt [M ⊙ pc − ] [Myr] [M ⊙ ] [M ⊙ ] [%] [%]8 26 . ± . . ± . × (3 . ± . × . ± . . ± . . ± . . ± . × (1 . ± . × . ± . . ± . Figure 3.
Cloud surface densities ( bottom row ), star formation rates ( second row ), virial parameters ( third row ), velocity dispersions( fourth row ), and radii ( top row ) for a set of 40 clouds. Half of the clouds were evolved with a low surface density characteristic of thebulk of the atomic ISM and the other half were evolved with a high surface density, characteristic of overdense regions in the ISM. Thetwo different choices of Σ res are marked at top. Accretion was shut off once 10 M ⊙ of material had been processed through the accretionflow. The ambient surface density, and thus the accretion history, strongly affects the resulting cloud evolution. ǫ issomewhat larger, approximately 8%. This can be entirelyattributed to the difference in lifetimes between the twosets of models. In the low surface density case, mostclouds are only able to survive one or two crossing timesand thus can only convert a small fraction of their massinto stars before they are dissociated or disrupted. Theclouds evolved with a high ambient surface density areable to survive for many crossing times and convert alarger fraction of their gas into stars. An even largerfraction is ejected via photoionization. However, for bothmodels, the star formation efficiency per free-fall time, ǫ ff = ˙ M ∗ M cl t ff (33)is low, around 2%. This is not surprising, as a low starformation efficiency per free-fall time is one of the basicassumptions of our model. Energetics of Star Formation Feedback VersusMass Accretion
GMCs exhibit highly supersonic turbulence. There isno agreement in the literature about what drives thesemotions, which numerical models of compressible MHDturbulence indicate should decay if left undriven. Someauthors suggest that the primary energy injection mech-anism is some sort of internal star formation feedbackprocess, such as protostellar outflows (Li & Nakamura2006; Wang et al. 2010), expanding H ii regions (Matzner2002), or supernovae (Mac Low & Klessen 2004). Oth-ers suggest that turbulence is driven externally via massinflows (Klessen & Hennebelle 2010). Comparing theamount of energy injected by different forms of star for-mation feedback, Fall et al. (2010) found that at typicalGMC column densities, the dominant stellar feedbackmechanism is H ii regions driven by the intense radia-tion fields emitted by massive star clusters. Using ourmodels, we can compare the importance of accretion rel-ative to H ii regions in the energy budget of GMCs.To find the total energy injected by accretion, we makeuse of our knowledge of the total energy of the cloud asa function of time. At the end of time step j we useequation (18) to calculate both the total cloud energy, E cl ,j , as well as what the cloud energy would have been ifwe had set ˙ M acc = 0 for that time step, E cl | ˙ M acc =0 . Thedifference, E acc , j = E cl ,j − E cl | ˙ M acc =0 (34)is the total energy added by accretion during that timestep. The total energy injected by accretion over thecloud’s lifetime is just the sum of the contributions ofeach time step, E acc = X j E acc ,j . (35)The energy injected by H ii region i , E H ii ,i , can befound by integrating the rate of energy injection by asingle H ii region with respect to time. This is, E H ii ,i = 1 . η E T ,i (cid:18) r m ,i R cl ,i (cid:19) / (36)where r m ,i is the radius of H ii region i when it mergeswith the parent cloud and R cl ,i is the radius of the cloud as a whole when H ii region i merged with the cloud.To find the total energy injected by H ii regions over thecloud’s lifetime, we simply sum up the contributions dueto individual H ii regions, E H ii = X i E H ii ,i . (37)The ratio |E H ii / E acc | indicates the relative importanceof star formation feedback to accretion driven turbulenceto the global energy budget of the cloud. If |E H ii / E acc | <
1, accretion dominates the energy injection; similarly if |E H ii / E acc | >
1, star formation feedback is the primarydriver of turbulence.The results of this comparison are plotted for bothchoices of the ambient surface density in Figure 4. Wefind that H ii regions and accretion contribute approxi-mately equal amounts of energy in the low surface den-sity runs, while accretion dominates in the high surfacedensity runs. In the low surface density runs, stochasticeffects can be important, particularly for clouds that donot last much longer than a crossing time. Thus, in someruns, star formation feedback can contribute significantlymore energy than accretion, while in others star forma-tion feedback is negligible. In the runs evolved with ahigh ambient surface density, star formation feedback issubdominant, although not completely negligible, in thevast majority of runs.It is worth pointing out that this result depends onthe precise value of ϕ we choose to evolve the cloudswith. If ϕ is lower, accretion contributes less energy, andstar formation can dominate the energy budget. If ϕ is higher, star formation becomes completely negligible,and the amount of kinetic energy injection is controlledby the mass accretion rate. Since clouds collapse when wechoose ϕ much lower than our fiducial value, and shocksin molecular gas tend to be strongly dissipative, we donot expect the ‘true’ value of ϕ to be much different thanour fiducial value. We thus conclude that one of threecases must hold. Star formation may be dominant, butonly marginally so. Accretion may also be dominant, butagain, only marginally. It is also possible that star for-mation and accretion contribute roughly equal amountsof energy. In all three cases, neither star formation oraccretion is truly negligible. OBSERVATIONAL COMPARISONS
Larson’s Laws
Giant molecular clouds are observed to obey threescaling relations, known as Larson’s laws (Larson 1981;Solomon et al. 1987; Bolatto et al. 2008). In their sim-plest form, Larson’s laws state: • The velocity dispersion scales with a power of thesize of the cloud. Subsequent observations haveshown that this power is about 0 . σ cl ∝ R . ). • The mass of the cloud scales with the square of theradius (constant Σ cl ). • Clouds are in approximate virial equilibrium ( α vir of order unity).These laws are not independent; any two imply theother. At a minimum, an acceptable theoretical modellobal Evolution of GMCs: Accretion 13 Figure 4.
The number of clouds plotted as a function of |E H ii / E acc | . In regions of low ambient surface density, accretion and starformation are in equipartition, while in regions of high ambient surface density, accretion dominates the energy budget. for GMCs should agree with both the scaling and thenormalization of the Larson scaling relations observed inreal clouds. We have already seen that clouds maintainapproximate virial equilibrium as well as roughly con-stant surface densities, but we have yet to see whetherthe normalization of the Larson scaling relations for ourmodels agrees with the observed Larson scaling relations. Equilibrium Surface Densities
GMCs, both in the Milky Way (Larson 1981;Solomon et al. 1987), and in nearby external galax-ies (Blitz et al. 2007; Bolatto et al. 2008) exhibit sur-prisingly little variation in surface density. For theSolomon et al. (1987) sample of Milky Way clouds, thiswas found to be h Σ cl i = 170 M ⊙ pc − . More recent andsensitive observations find lower values, closer to h Σ cl i =50 M ⊙ pc − in the Milky Way (Heyer et al. 2009) andin the Large Magellanic Cloud (Hughes et al. 2010), al-though these latter estimates depend on a highly uncer-tain correction for non-LTE line excitation and the COto H conversion factor, respectively. Using heterogenousdata from several nearby galaxies, Bolatto et al. (2008)attempted to extract cloud properties in a uniform man-ner and found a typical surface density of 85 M ⊙ pc − but with significant variation from galaxy to galaxy.Variations in the mean GMC surface density are seenwhen comparing samples from different galaxies. How-ever, within a single galaxy there is little variation(Blitz et al. 2007). These variations are usually at-tributed to differences in the CO-to-H conversion factorfrom galaxy to galaxy (Bolatto et al. 2008), a quantitywhich may depend on metallicity and the interstellar ra-diation field (Glover et al. 2010) as well as variations inturbulent pressure and radiation field in the ambient in-terstellar medium. In our runs, we also recover roughlyconstant surface densities (see the second row from the bottom of Figure 3).In Figure 5, we have reproduced a figure fromBlitz et al. (2007) that depicts observational results forCO luminosities and cloud radii for a sample of cloudsin the outer Milky Way as well as from several samplesof extragalactic GMCs. To compare against this com-pendium of results, we calculate CO luminosities for ourmodel clouds by assuming a constant CO-to-H conver-sion factor, L CO = M cl . ⊙ K km s − pc (38)as in Rosolowsky & Leroy (2006). This formula ac-counts for the presence of helium and assumes aconstant H to CO conversion factor, X CO = 4 × cm − ( K km s − ) − , twice the value derived formolecular clouds within the Solar circle using observa-tions of gamma-ray emission (Strong & Mattox 1996;Abdo et al. 2010). We choose this value to be consis-tent with Blitz et al. (2007), who find, using this valueof X CO , that all of the GMCs in their sample have virialmasses comparable to the masses implied by their COluminosity to within a factor of two.With our fiducial initial conditions, model clouds in oursample begin their lives in the bottom left hand cornerof Figure 5, at R cl ≈
10 pc. As they accrete and ex-pand, clouds move towards the upper right hand corner.Clouds end their evolution either through disruption bya single H ii region or by passing below the moleculardissociation threshold, indicated by a dashed line in Fig-ure 5. Offsets in the distribution of column densities fromgalaxy to galaxy and from the simulated clouds can beattributed to variations in X CO and uncertainty in iden-tifying a unique radius for observed clouds (Blitz et al.2007) that have nonzero obliquity (Bertoldi & McKee1992). Accounting for variations in X CO , there is striking4 Goldbaum et al. Figure 5.
Cloud CO luminosity plotted as a function of cloudradius. CO luminosities are found by assuming X CO = 4 × cm − ( K km s − ) − . Solid lines of constant surface densityare plotted for 10, 100, and 1000 M ⊙ pc − for reference. Thedashed line of constant surface density corresponds to our assumeddissociation threshold. The outputs from a set of 2000 runs wereused, with Σ res = 8 and 16 M ⊙ pc − and M res = 6 × M ⊙ .Colors indicate the amount of time model clouds tend to occupya position in L CO – R cl parameter space. Symbols denote ob-served CO luminosities and cloud radii for galactic ( points ) andextragalactic ( open shapes ) GMCs. See Blitz et al. (2007) and ref-erences therein for details of the observations. agreement between the observed distribution of molecu-lar clouds and our sample of simulated clouds.The models exhibit a kink in their evolution whenthe reservoir is exhausted and accretion is shut off.For this reason, there are no clouds with L CO > . K km s − pc . Once accretion is shut off, theclouds decrease in mass for the remainder of their evo-lution. This kink is somewhat artificial since we haveassumed a fixed reservoir mass and a smooth accretionhistory. A more sophisticated model for the reservoir in-cluding a range of reservoir masses would exhibit a con-tinuous spectrum of kinks, broadening the region of pa-rameter space explored by the models, particularly for L CO < . K km s − pc .The models also exhibit two distinct favored strips ofparameter space along which they tend to evolve. Thiscorresponds to the two different equilibrium column den-sities picked out by the two different choices of Σ res . Thisbehavior is clearly seen in the second panel from the bot-tom of Figure 3. The fact that Σ cl is sensitive to Σ res follows from dimensional analysis. Linewidth-Size Relation
We next compare our simulated clouds with thelinewidth-size relation observed to hold among GMCs asa population (Bolatto et al. 2008). We are able to repro-duce the power law, scatter, and the rough normalizationin the observed linewidth-size relation. This conclusionis unsurprising, since we have already seen that our simu-lated clouds maintain roughly constant virial parameters
Figure 6.
Cloud velocity dispersion plotted as a function of cloudradius. The dash-dotted line is the galactic linewidth-size relationfound by Solomon et al. (1987), σ v = 0 . R . . Symbols and colorcoding are the same as in Figure 5. and surface densities as they evolve. It is worth notingthat for our simplified model for the environment of aGMC, that the linewidth-size relation corresponds to anage sequence. Clouds that live towards the left-hand sideof the diagram are younger than clouds that live towardsthe right. It is possible that this conclusion is an arti-fact of choosing a single reservoir mass. Clouds accretingfrom a population of reservoirs with a continuous spec-trum of masses may blur this effect somewhat. We planto revisit this in future work in which we will model theglobal ISM of a galaxy simultaneously with the evolutionof a population of GMCs.There is a small offset when comparing the locus ofextragalactic and outer Milky Way clouds with our mod-els, although there is good agreement between our mod-els and the scaling found by Solomon et al. (1987). Fora subset of the observational sample, particularly theSMC clouds, it is possible that the metallicity of thegas in the clouds is so low that CO is no longer agood tracer of the bulk of the molecular gas (Leroy et al.2007). Since our models assume perfect sphericity, andthe observed radius of a prolate or oblate spheroid willalways be smaller than the corresponding spherical ra-dius (Bertoldi & McKee 1992), it is also possible thatthe radii predicted by our models overpredict the corre-sponding observed cloud radius by 0.1 or 0.2 dex. Lastly,it could be the that we overpredict the various pres-sures due to photoionization and accretion by assumingspherical symmetry. In reality, the wind and accretionram pressure may not necessarily be perfectly sphericallysymmetric, leading to a reduction in the overall confiningpressure and an increase in the radius. Evolutionary Classification
The Large Magellanic Cloud is home to one of the best-studied samples of GMCs in any galaxy. The LMC’slobal Evolution of GMCs: Accretion 15disk-like geometry and face-on orientation offers littleambiguity in distance measurements, with the most accu-rate measurements giving d LMC = 50 . CO ( J = 1 →
0) surveysand high-resolution followup from the MAGMA CO( J = 1 →
0) survey (Hughes et al. 2010) have mappedthe molecular content of the entire disk of the LMC andidentified 272 clouds that together contain 5 × M ⊙ of molecular gas. When combined with multiwavelengtharchival observations of star formation indicators, theseCO data constitute a snapshot in the evolution and starformation history of a population of GMCs.Kawamura et al. (2009) used the NANTEN CO J =(1 →
0) data, along with complementary H α photom-etry (Kennicutt & Hodge 1986), radio continuum mapsat 1 .
4, 4 .
8, and 8 . <
10 Myr) clusters extractedfrom
U BV photometry (Bica et al. 1996) to investigatethe ongoing star formation within GMCs in the LMC.These authors found a strong tendency for H ii regionsand young clusters to be spatially correlated with GMCs.Using this association, the GMCs in their sample wereseparated into three types. Type 1 GMCs are defined tobe starless in the sense that they are not associated withdetectable H ii regions or young clusters, Type 2 GMCsare associated with H ii regions, but not young clustersin the cluster catalog, and Type 3 GMCs are associatedwith both H ii regions and young clusters. 24% of theNANTEN sample were classified as Type 1, 50% as Type2, and 26% as Type 3.Assuming that GMCs and clusters are formed in steadystate and assuming that young clusters not associatedwith GMCs are associated with GMCs that have dis-sipated, one can infer from the NANTEN populationstatistics that GMCs spend 6 Myr in the Type 1 phase,13 Myr in the Type 2 phase, 7 Myr in the Type 3 phase,and then dissipate within 3 Myr. This accounting im-plies GMC lifetimes of approximately 20to30 Myr. Insupport of the claim that the GMC classification schemeconstitutes an evolutionary sequence, the authors notethat among the resolved GMCs in the NANTEN survey,Type 3 GMCs are on average more massive, have largerturbulent line widths, and have larger radii. However,there is significant scatter in the Type 3 GMC sampleand the mass and size evolution are well within theirerror bars.In order to correct for extinction, which might obscureH α emitting H ii regions, radio continuum maps at three,well-separated frequencies were used to identify obscuredH ii regions via their flat spectral slopes. However, noH ii regions were identified in the radio continuum datathat were not present in the H α maps, leading the au-thors to conclude that the H α data was unaffected byobscuration. No similar analysis was performed to esti-mate obscuration of young star clusters. No attempt wasmade to correct for the varying sensitivities in the dif-ferent radio maps, allowing for the possibility that someH ii regions were detected at 1.4 Ghz but below the sen-sitivity limit at 4.8 and 8.6 Ghz.There are several observational biases inherent in theGMC classification scheme described above. The first is the probable existence of star clusters and H ii regionslocated either behind or within giant molecular cloudsfrom our viewpoint. High dust extinction along thesesightlines would mask some young clusters from detec-tion in the Bica et al. (1996) star cluster sample. Thiscould lead to an overestimate of Type 2 GMCs relativeto Type 3 GMCs. Another possible bias is the use ofthe Bica et al. (1996) star cluster catalog. Clusters inthis catalog were targeted for U BV photometry basedon brightness and association with emission nebulae. Itis possible that some young clusters were missed in thiscatalog and no attempt is made by Kawamura et al. tocorrect for the completeness of the cluster catalog. Thiswould also lead to an overestimate of Type 2 GMCs rel-ative to Type 3 GMCs.In order to make a quantitative comparison be-tween our models and the evolutionary classification ofKawamura et al. (2009), we employ a few simple pre-scriptions to generate synthetic V band, H α , and ra-dio continuum photometry for the clusters and H ii re-gions in our simulated clouds. First, we calculate the V -band luminosity of our simulated clusters using thesynthetic photometry of Lejeune & Schaerer (2001). For5 M ⊙ ≤ M ∗ ≤
15 M ⊙ , the photometry is based on theevolution tracks of Schaerer et al. (1993). For massivestars, M ∗ ≥
15 M ⊙ , the synthetic photometry is basedon the high mass-loss models of Meynet et al. (1994). Forboth sets of synthetic photometry, we assume Z = 0 . L V for these stars by only considering the mainsequence evolution and taking L V = h L V i MS , the meanluminosity on the main sequence. Since the stellar evo-lutionary tracks give the luminosity at a discrete set ofmasses, we interpolate by fitting a broken power law be-tween stellar masses with evolutionary tracks.We calculate the H α luminosity of our model H ii re-gions via (McKee & Williams 1997), L H α = 1 . × S erg s − , (39)where S is the ionizing luminosity of the central starcluster in units of 10 photon s − . This is larger thanthe empirical relation by a factor of 1 .
37 to correct for theabsorption of ionizing radiation by dust grains. Lastly,we find the radio continuum luminosity of our simulatedH ii regions via (Condon 1992), L ν = 1 . × S (cid:18) ν (cid:19) erg s − Hz − (40)We also account for the reduction in flux from GMCs thatwould be larger than the beam size used by Dickel et al.(2005) by performing a geometric correction and assum-ing that all H ii regions are placed at the center of themodel beam.To assign our simulated GMCs an evolutionary classi-fication, we extract the ionizing luminosity and V bandluminosity of the brightest cluster in the GMC as a func-tion of time. Using the ionizing flux, we calculate theexpected H α and radio continuum luminosity via equa-tions (39) and (40) respectively. For the V band and Hα luminosity, we correct for the foreground extinction, A V = 0 .
25, towards the LMC (Schlegel et al. 1998), ig-noring extinction internal to the LMC but external tothe cloud under consideration. To identify the GMC as6 Goldbaum et al.
Figure 7.
Fraction of GMC lifetime spent as a Type 1, 2, and 3GMC as defined by Kawamura et al. (2009) ( diamonds, top row ),and average age of GMCs in each classification bin ( bottom row ).The asterisks indicate the median of the GMC type probabilitydistribution functions generated using a Monte Carlo analysis de-scribed in the text. The error bars encompass the 10th to 90thpercentile interval of the probability distribution functions. either Type 2 or Type 3, we use either an H α luminos-ity cutoff or a radio continuum flux cutoff correspond-ing to the detection limits quoted by (Kawamura et al.2009). For H α , we say that an H ii region is detected if L Hα ≥ erg s − . For radio continuum, we say an H ii region is detected if the radio flux at a distance of 50 kpcwould be greater than 0 . − for a beam size of20 ′′ at 4 . α or radio continuum as well as if L V ≥ . × L ⊙ ,the completeness limit quoted by Bica et al. (1996) forthe young cluster sample.Since there are bound to be clusters that are locatedboth behind and within clouds along our line of sight, wealso correct the synthetic V -band and H α photometryfor extinction by the GMC. This is done by assumingthat star clusters form at random locations within clouds.We calculate the optical depth to the star cluster via τ ν = κ ν ξ µ R cl where κ ν is the mean dust opacity throughthe cloud, ξ µ is the depth into the cloud in units of thecloud radius and µ is the angle between the normal to thesurface of the cloud and the line of sight (Krumholz et al.2008). For this purpose, we use a Milky Way extinctioncurve and assume an LMC dust-to-gas ratio whereby acolumn of 1 g cm − corresponds to A V = 107. In thevisual passbands we are concerned with here, the MilkyWay and LMC extinction curves are nearly identical.We have run a set of 2000 cloud models, 1000 eachevolved using two different values of the ambient surfacedensity, Σ res = 8 and 16 M ⊙ pc − . All other parametersare as in Table 1. The resulting cloud models encompassthe entire observed range in cloud masses reported byFukui et al. (2008).To directly compare to the observed distribution ofGMC types we perform simulated observations using aMonte Carlo scheme. Since the observations are inher-ently weighted by the GMC mass function, we first gen-erate cloud masses by drawing from a powerlaw GMCmass spectrum with a slope of -1.6, a minimum massof 5 × M ⊙ and a maximum mass of 5 × M ⊙ (Fukui & Kawamura 2010). Once a mass is generated, Figure 8.
GMC classification as a function of time for a selectionof model clouds. we find all time steps where model clouds have masseswithin 0.1 dex of the randomly selected mass. Withinthis sample of time steps, we calculate f , f , and f ,the fraction of Type 1, Type 2, and Type 3 GMCs, re-spectively. At the same time, we calculate t , t , and t the average age of clouds in each GMC Type bin. Wegenerate 10 Monte Carlo realizations, from which weconstruct probability distributions for f , f , f , t , t ,and t .The results of this comparison are presented in the toprow of Figure 7. In the figure, the lines connect the me-dian of the Monte Carlo probability distributions whilethe error bars encompass the 10th and 90th percentile.We are able to reproduce the observed distribution ofType 1, 2 and 3 GMCs as observed by Kawamura et al.(2009). In particular, using both detection limits, we findthat the majority of clouds are detected as Type 2 GMCs,and relatively fewer clouds are detected as Type 1 and 3GMCs. Interestingly, in the bottom panel of the figure,we find that, on average, the GMC classification schemedoes constitute an age sequence in that Type 2 GMCstend to be somewhat older than Type 1 clouds. Type 3GMCs in turn tend to be older than Type 2 clouds. Onthe other hand, the spread in cloud ages within each binis well within the error bars, indicating that the GMCtype classification is not necessarily a strict evolutionarysequence: older clouds can be classified as Type 1 andyounger clouds can be classified as Type 3.This can be seen more clearly in Figure 8, where weplot the classification of a selection of GMCs as a func-tion of time. We see that for some clouds, the classifica-tion scheme does represent an evolutionary sequence. Inthese runs, thw cloud starts as a Type 1 GMC, beginsforming star clusters, evolves into a Type 2 GMC, andthen forms massive OB association and becomes a Type3 GMC. However, we also see that a cloud can quicklyform a massive OB association and be classified as a Type3 GMC early on and only later be classified as a Type2 GMC. Alternatively, a cloud may happen to not formany massive clusters late in its evolution, causing a mas-sive and old cloud to be identified as a Type 2 or 1 GMClate in its evolution. Finally, there are clouds which ex-hibit no discernible pattern in their histories, more orless randomly transitioning between GMC classificationsthroughout their lives. This may explain the presence ofmassive ( ∼ M ⊙ ) Type 1 “young” GMCs in the sam-ples of Kawamura et al. (2009) and Hughes et al. (2010). CAVEATS AND LIMITATIONS lobal Evolution of GMCs: Accretion 17
Implications of the Assumption of Homology
Assuming that clouds evolve homologously is the mainlimitation of our model. We make the assumption of ho-mology to significantly simplify the equations governingthe evolution of the cloud. Given the assumption of ho-mology, our equations of motion follow rigorously fromthe local equation of momentum and energy conserva-tion. A more complex cloud structure destroys the rela-tive simplicity of the model and would require computa-tionally expensive hydrodynamical simulations to modelaccretion in detail.Homology constrains the cloud to always maintain thesame shape and degree of central concentration. This isequivalent to setting time derivates of the structure con-stant, a I , to zero. This might be a problem if changes ofthe moment of inertia of clouds occur primarily by chang-ing the shape of the cloud rather than through overall ex-pansion or contraction. It might also be a problem if theaccretion flow does not in reality differentially mix withthe cloud. If the accretion flow is anisotropic or if it can-not fall to the central regions of the cloud before it mixes,the cloud may become less centrally condensed and wemay overestimate the kinetic energy injection since ma-terial cannot fall to the bottom of the cloud potentialwell.While we cannot resolve the dynamical effect ofchanges in the shape of the cloud, we can resolve changesin the size of the cloud. Given the observed Larsonscaling relations, we expect more massive clouds to belarger, implying that clouds must expand as they accretemass. If clouds do form via gravitational instability, theymust accrete significant mass. We do resolve this behav-ior in our models. We caution that our virial analysismost readily describes a relaxed system. We may be do-ing a poor job of resolving phenomena that occur on atimescale comparable to the crossing time. Magnetic fields in the Atomic Envelope
One key assumption we made in deriving the global en-ergy equation was that the atomic envelope contributesnegligibly to the total magnetic energy associated withthe cloud. That is, we did not include the magnetic fieldwhen calculating E amb in equation (C6). The motiva-tion for this assumption is based on comparisons of ob-servations of magnetic fields in dense molecular clumps,where it is possible to measure the magnetic field di-rectly via the Zeeman effect in lines of OH and CN(Troland & Crutcher 2008; Falgarone et al. 2008), to themagnetic fields in the atomic ISM, measured via theZeeman effect in the 21 cm line of neutral Hydrogen(Heiles & Troland 2005). These studies consistently findthat the magnetic field strength is significantly elevatedin the dense molecular gas. However, most of the volumeof a GMC is occupied by diffuse molecular gas with lowOH abundance — frustrating efforts to measure magneticfields via observations of OH in emission. Thus, withoutdirect observations of magnetic fields in the diffuse molec-ular gas, we cannot know whether the magnetic fieldsin the atomic envelope are weak or strong compared tothe magnetic field strength in the bulk of a GMC. Whilepreliminary observations by Troland et al. (2011, privatecommunication) indicate that magnetic fields in the dif-fuse molecular component are somewhat stronger than in the atomic envelope, this question has yet to be settled.It is possible that we underestimate the net magneticenergy due to the presence of the cloud.Another possible problem with our treatment of mag-netic fields is that we assume the mass-to-flux ratio re-mains constant throughout the evolution of the cloud.This is equivalent to assuming a fixed value for η B .This assumption might be invalid if accreted mate-rial flows preferentially along magnetic field lines orif ambipolar diffusion can act act on timescales com-parable to the cloud dynamical time. While mea-surements (Li et al. 2006b) and theoretical esimates(McKee & Ostriker 2007) of the mass-to-flux ratio ofmolecular clouds find that clouds should be marginallysupercritical, with mass-to-flux ratios of order unity, thetime-evolution and cloud-to-cloud variation in mass-to-flux ratios are poorly constrained. Validity of Larson’s Laws
It has been suggested that the observed con-stancy of GMC surface densities is an artifact(see e.g. Kegel 1989; Vazquez-Semadeni et al.1997; Ballesteros-Paredes & Mac Low 2002;Ballesteros-Paredes 2006; Ballesteros-Paredes et al.2011). Since the physical structure we assume for theclouds in our model implicitly assumes that the Larsonrelations are valid, this may imply that our models donot correspond well to real GMCs.The argument that the Larson laws are a product of aselection effect usually proceeds as follows. If one looksfor overdensities of a particular size in a simulation of theturbulent ISM, one will find that the selected clouds havea wide distribution of masses, implying a large spreadin cloud-to-cloud surface density. Since observers detectclouds using CO as a tracer and at low surface densitythe CO abundance is not high enough to be detectablein emission, observers will never find low surface densityclouds. This implies that at a fixed radius, real cloudsshould have more variation in mass than a naive inter-pretation of the CO observations would suggest.This argument misses two key aspects of the observedproperties of GMCs. The first is that it cannot explainthe lack of clouds at high gas surface densities. For thesame reason that we should not be able to see diffuseclouds, we should very easily be able to see compact, highsurface density CO clouds. The fact that these cloudsdon’t exist implies something important about molecu-lar cloud structure. The second argument is that the lackof low surface density clouds does not imply that molec-ular clouds can exist at all surface densities but merelythat the molecular clouds dissociate once they becomeoptically thin to the ambient ultraviolet radiation field.While diffuse atomic clouds certainly exist, these cloudsdo not form stars (Krumholz et al. 2011).Since high surface density GMCs are not observed inthe local universe and low surface density clouds arenot molecular and thus not GMCs by definition, the ob-served lack of variation in GMC surface densities mustbe a real property of the clouds. This has been con-firmed with very detailed dust extinction measurementsof nearby star forming clouds, where an exquisitely tightmass-radius relation is observed (Lombardi et al. 2010)and in extragalactic studies where little variation is seenwhen comparing the mass-radius relation from galaxy8 Goldbaum et al.to galaxy (Bolatto et al. 2008). Taken together, the ev-idence seems to imply that the Larson relations are aproperty of the structure of GMCs and are not due to aselection effect. CONCLUSIONS
In this paper, we have presented semianalytic dynami-cal models for the evolution of giant molecular clouds un-dergoing both mass accretion and star formation. Thesemodels are able to capture the evolution of individualGMCs from their growth and the onset of massive starformation, until their dispersal via an energetic H ii re-gion or through the combined action of accretion andstar formation. We are able for the first time to syn-thesize galactic populations of GMCs whose propertiescorrespond closely to the observed properties of GMCsin the Milky Way and nearby external galaxies. We haveshown that clouds in low surface density environmentsgenerally disperse within a few crossing times, beforethey can accrete all of the gas in their reservoir. Atthe same time, clouds in high surface density environ-ments do accrete all of the gas in their reservoirs andtend to be larger and more massive. We have also shownthat mass accretion can contribute a significant frac-tion of the total energy available for turbulent driving.Lastly, we generate synthetic cluster observations andcompare against the evolutionary classification schemeof Kawamura et al. (2009), finding good agreement whenwe correct for selection effects and systematic biases in-herent in the observations. We conclude that, on aver-age, the evolutionary classification scheme correspondsto an age sequence but is not a good predictor for theevolutionary state of isolated clouds.We would like to thank Leo Blitz and Adam Leroyfor making available the Milky Way data used in Fig-ures 5 and 6. We would also like the thank an anony-mous referee for their helpful comments which spurreda significant improvement in the quality of the paper.NJG is supported by a Graduate Research Fellowshipfrom the National Science Foundation. MRK acknowl-edges support from: an Alfred P. Sloan Fellowship; NSFgrants AST-0807739 and CAREER-0955300; and NASAthrough Astrophysics Theory and Fundamental Physicsgrant NNX09AK31G and a Spitzer Space Telescope The-oretical Research Program grant. CDM’s research issupported by an NSERC discovery grant and an On-tario Early Researcher Award. The research of CFM issupported in part by NSF grant AST-0908553 and by aSpitzer Space Telescope Theoretical Research Programgrant. REFERENCESAbel, T., Bryan, G. L., & Norman, M. L. 2002, Science, 295, 93Abdo, A. A., Ackerman, M., Ajello, M., et al. 2010, ApJ, 710, 133Alves, D. R. 2004, New Astronomy Review, 48, 659Ballesteros-Paredes, J., Hartmann, L., & V´azquez-Semadeni, E.1999, ApJ, 527, 285Ballesteros-Paredes, J., & Mac Low, M.-M. 2002, ApJ, 570, 734Ballesteros-Paredes, J. 2006, MNRAS, 372, 443Ballesteros-Paredes, J., Hartmann, L. W., V´azquez-Semadeni, E.,Heitsch, F., & Zamora-Avil´es, M. A. 2011, MNRAS, 411, 65Bertoldi, F., & McKee, C. F. 1992, ApJ, 395, 140 Bica, E., Claria, J. J., Dottori, H., Santos, J. F. C., Jr., & Piatti,A. E. 1996, ApJS, 102, 57Blitz, L., Fukui, Y., Kawamura, A., Leroy, A., et al. 2007, inProtostars and Planets V, ed. B. Reipurth, D. Jewitt, & K.Keil (Tucson, AZ; Univ. Arizona Press), 81Blitz, L., & Shu, F. H. 1980, ApJ, 238, 148Bigiel, F., Leroy, A., Walter, F., et al. 2008, AJ, 136, 2846Bolatto, A. D., Leroy, A. K., Rosolowsky, E., Walter, F., & Blitz,L. 2008, ApJ, 686, 948Brunt, C. M., & Heyer, M. H. 2002, ApJ, 566, 276Brunt, C. M., Heyer, M. H., & Mac Low, M.-M. 2009, A&A, 504,883Condon, J. J. 1992, ARA&A, 30, 575Dickel, J. R., McIntyre, V. J., Gruendl, R. A., & Milne, D. K.2005, AJ, 129, 790Dobbs, C. L., Burkert, A., & Pringle, J. E. 2011, MNRAS, 413,2935Elmegreen, B. G., & Scalo, J. 2004, ARA&A, 42, 211Engargiola, G., Plambeck, R. L., Rosolowsky, E., & Blitz, L.2003, ApJS, 149, 343Evans, N. J., et al. 2009, ApJS, 181, 321van Dishoeck, E. F., & Black, J. H. 1988, ApJ, 334, 771Falgarone, E., Troland, T. H., Crutcher, R. M., & Paubert, G.2008, A&A, 487, 247Fall, S. M., Krumholz, M. R., & Matzner, C. D. 2010, ApJ, 710,L142Fukui, Y., Kawamura, A., Wong, T., et al. 2008, ApJS, 178, 56Fukui, Y., Kawamura, A., Minamidani, T., et al. 2009, ApJ, 705,144Fukui, Y., & Kawamura, A. 2010, ARA&A, 48, 547Glover, S. C. O., Federrath, C., Mac Low, M.-M., & Klessen,R. S. 2010, MNRAS, 404, 2Goldreich, P., & Kwan, J. 1974, ApJ, 189, 441Heiles, C., & Troland, T. H. 2005, ApJ, 624, 773Heyer, M. H., Carpenter, J. M., & Snell, R. L. 2001, ApJ, 551, 852Heyer, M. H., & Brunt, C. M. 2004, ApJ, 615, L45Heyer, M., Krawczyk, C., Duval, J., & Jackson, J. M. 2009, ApJ,699, 1092Hughes, A., Staveley-Smith, L., Kim, S., Wolleben, M., &Filipovi´c, M. 2007, MNRAS, 382, 543Hughes, A., Wong, T., Ott, J., et al. 2010, MNRAS, 406, 2065Hunter, C. 1977, ApJ, 218, 834Imara, N., Bigiel, F., & Blitz, L. 2011, ApJ, 732, 79Kalberla, P. M. W., & Dedes, L. 2008, A&A, 487, 951Kawamura, A., Mizuno, Y., Minamidani, T., et al. 2009, ApJS,184, 1Kegel, W. H. 1989, A&A, 225, 517Kennicutt, R. C., Jr., & Hodge, P. W. 1986, ApJ, 306, 130Kennicutt, R. C., Jr., Calzetti, D., Water, F., et al. 2007, ApJ,671, 333Kennicutt, R. C., Jr. 1998, ApJ, 498, 541Kim, W.-T., Ostriker, E. C., & Stone, J. M. 2002, ApJ, 581, 1080Kim, W.-T., Ostriker, E. C., & Stone, J. M. 2003, ApJ, 599, 1157Kim, W.-T., & Ostriker, E. C. 2006, ApJ, 646, 213Klessen, R. S., & Hennebelle, P. 2010, A&A, 520, A17Kroupa, P. 2002, Science, 295, 82Krumholz, M. R., & McKee, C. F. 2005, ApJ, 630, 250Krumholz, M. R., Matzner, C. D., & McKee, C. F. 2006, ApJ,653, 361Krumholz, M. R., & Tan, J. C. 2007, ApJ, 654, 304Krumholz, M. R., McKee, C. F., & Tumlinson, J. 2008, ApJ, 689,865Krumholz, M. R., & Matzner, C. D. 2009, ApJ, 703, 1352Krumholz, M. R., Leroy, A. K., & McKee, C. F. 2011, ApJ, 731,25Lada, C. J., Lombardi, M., & Alves, J. F. 2010, ApJ, 724, 687Larson, R. B. 1981, MNRAS, 194, 809Lejeune, T., & Schaerer, D. 2001, A&A, 366, 538Leroy, A., Bolatto, A., Stanimirovic, S., et al. 2007, ApJ, 658,1027Leroy, A. K., Walter, F., Brinks, et al. 2008, AJ, 136, 2782Li, Y., Mac Low, M.-M., & Klessen, R. S. 2006, ApJ, 639, 879Li, H., Griffin, G. S., Krejny, M., et al. 2006, ApJ, 648, 340Li, Z.-Y., & Nakamura, F. 2006, ApJ, 640, L187Li, Z.-Y., Wang, P., Abel, T., & Nakamura, F. 2010, ApJ, 720,L26Lombardi, M., Alves, J., & Lada, C. J. 2010, A&A, 519, L7 lobal Evolution of GMCs: Accretion 19
Lopez, L. A., Krumholz, M. R., Bolatto, A. D., Prochaska, J. X.,& Ramirez-Ruiz, E. 2011, ApJ, 731, 91Mac Low, M.-M., Klessen, R. S., Burkert, A., & Smith, M. D.1998, Physical Review Letters, 80, 2754Mac Low, M.-M., & Klessen, R. S. 2004, Reviews of ModernPhysics, 76, 125Matzner, C. D., & McKee, C. F. 2000, ApJ, 545, 364Matzner, C. D. 2001, in ASP Conf. Proc. 243, From Darkness toLight: Origin and Evolution of Young Stellar Clusters, ed. T.Montmerle, P. Andr´e (San Francisco, CA: ASP), 757Matzner, C. D. 2002, ApJ, 566, 302McKee, C. F. 1989, ApJ, 345, 782McKee, C. F., Zweibel, E. G., Goodman, A. A., & Heiles, C.1993, in Protostars and Planets III, ed. E. H. Levy, & J. I.Lunine (Tucson, AZ: Univ. Arizona Press), 327McKee, C. F. 1999, in NATO ASIC Proc. 540: The Origin ofStars and Planetary Systems, ed. C. J. Lada & N. D. Kylafis(Dordrecht: Kluwer), 29McKee, C. F., & Ostriker, E. C. 2007, ARA&A, 45, 565McKee, C. F., & Tan, J. C. 2003, ApJ, 585, 850McKee, C. F., & Williams, J. P. 1997, ApJ, 476, 144McKee, C. F., & Zweibel, E. G. 1992, ApJ, 399, 551Meynet, G., Maeder, A., Schaller, G., Schaerer, D., &Charbonnel, C. 1994, A&AS, 103, 97Mizuno, N., Rubio, M., Mizuno, A., et al. 2001, PASJ, 53, L45Norman, C., & Silk, J. 1980, ApJ, 238, 158Murray, N. 2011, ApJ, 729, 133Onodera, S., et al. 2010, ApJ, 722, L127Ossenkopf, V., & Mac Low, M.-M. 2002, A&A, 390, 307Parravano, A., Hollenbach, D. J., & McKee, C. F. 2003, ApJ, 584,797Roman-Duval, J., Jackson, J. M., Heyer, M., Rathborne, J., &Simon, R. 2010, ApJ, 723, 492Rosolowsky, E. 2005, PASP, 117, 1403Rosolowsky, E. 2007, ApJ, 654, 240Rosolowsky, E., & Leroy, A. 2006, PASP, 118, 590Schaerer, D., Meynet, G., Maeder, A., & Schaller, G. 1993,A&AS, 98, 523Schruba, A., Leroy, A. K., Walter, F., Sandstrom, K., &Rosolowsky, E. 2010, ApJ, 722, 1699Schruba, A., Leroy, A. K., Water, F., et al. 2011, arXiv:1105.4605 Shu, F. H. 1977, ApJ, 214, 488Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500,525Solomon, P. M., Rivolo, A. R., Barrett, J., & Yahil, A. 1987, ApJ,319, 730Stahler, S. W., Shu, F. H., & Taam, R. E. 1980, ApJ, 241, 637Stark, A. A., & Lee, Y. 2006, ApJ, 641, L113Stone, J. M., Ostriker, E. C., & Gammie, C. F. 1998, ApJ, 508,L99Strong, A. W., & Mattox, J. R. 1996, A&A, 308, L21Tamburro, D., Rix, H.-W., Walter, et al. 2008, AJ, 136, 2872Tan, J. C., & McKee, C. F. 2004, ApJ, 603, 383Tan, J. C., Krumholz, M. R., & McKee, C. F. 2006, ApJ, 641,L121Tasker, E. J., & Tan, J. C. 2009, ApJ, 700, 358Tasker, E. J. 2011, ApJ, 730, 11Thilker, D. A., Braun, R., & Walterbos, R. A. M. 2002, in ASPConf. Ser. 276, Seeing Through the Dust: The Detection of HIand the Exploration of the ISM in Galaxies, ed, A. R. Taylor,T. L. Landecher, & A. G. Willis (San Francisco, CA: ASP), 370Troland, T. H., & Crutcher, R. M. 2008, ApJ, 680, 457van Dishoeck, E. F., & Black, J. H. 1988, ApJ, 334, 771V´azquez-Semadeni, E., G´omez, G. C., Jappsen, A. K., et al. 2007,ApJ, 657, 870Vazquez-Semadeni, E., Ballesteros-Paredes, J., & Rodriguez,L. F. 1997, ApJ, 474, 292V´azquez-Semadeni, E., Col´ın, P., G´omez, G. C.,Ballesteros-Paredes, J., & Watson, A. W. 2010, ApJ, 715, 1302Wang, P., & Abel, T. 2008, ApJ, 672, 752Wang, P., Li, Z.-Y., Abel, T., & Nakamura, F. 2010, ApJ, 709, 27Williams, J. P., & McKee, C. F. 1997, ApJ, 476, 166Wise, J. H., & Abel, T. 2007, ApJ, 665, 899Wolfire, M. G., Hollenbach, D., & McKee, C. F. 2010,arXiv:1004.5401Wong, T., & Blitz, L. 2002, ApJ, 569, 157Wong, T., et al. 2009, ApJ, 696, 370Yang, C.-C., Gruendl, R. A., Chu, Y.-H., Mac Low, M.-M., &Fukui, Y. 2007, ApJ, 671, 374Zuckerman, B., & Evans, N. J., II 1974, ApJ, 192, L149APPENDIX A. DERIVATION OF THE VIRIAL THEOREM FOR AN ACCRETING CLOUD
Here we derive an equation governing the virial balance of a cloud that is simultaneously forming stars and accretingmaterial. This is a generalization of the analysis of Paper I and McKee & Zweibel (1992). We refer the reader to thosepapers for details that are unrelated to the accretion flow and to § V vir with bounding surface S vir . We assumethat V vir is sufficiently large to contain the cloud at all times. Material within the virial volume is apportioned intothree components: virial material, a gaseous reservoir, and material in a photoionized wind. Locally, each componentsatisfies its own continuity equation, ∂ρ∂t = − ∇ · ρ v + ˙ ρ (A1) ∂ρ res ∂t = − ∇ · ρ res v res − ˙ ρ acc (A2) ∂ρ w ∂t = − ∇ · ρ w v w − ˙ ρ ej . (A3)These equations are coupled via the source and sink terms,˙ ρ = ˙ ρ acc + ˙ ρ ej . (A4)We assume that accretion can only transport material onto the cloud and that the wind can only carry material awayfrom the cloud, implying ˙ ρ acc ≥ ρ ej ≤ ∂∂t ( ρ v ) = − ∇ · (Π − T M ) + ρ g + F ∗ + ˙ ρ ej ( v + v ′ ej ) + ˙ ρ acc v res (A5)where Π = P th I + ρ vv − π is the gas pressure tensor, π is the viscous stress tensor, T M = [ BB − (1 / B I ] / (4 π ) is0 Goldbaum et al.the Maxwell stress tensor, B is the magnetic field, I is the unit tensor, g is the gravitational force per unit mass, and F ∗ + ˙ ρ ej ( v + v ′ ej ) is the local body force due to the interaction between expanding H ii regions and the cloud. Wewrite the stellar forcing term in this form so as to separate the random component, F ∗ , from the spherically symmetriccomponent ˙ ρ ej ( v + v ′ ej ). F ∗ is a function of time and position in the cloud and depends on the precise location and history of massive starformation. Due to the nature of supersonic turbulence, we expect turbulent overdensities to be randomly scatteredthroughout the cloud. Thus, we expect F ∗ to be randomly oriented with respect to the position vector, implying R V vir r · F ∗ dV = 0. While F ∗ is randomly oriented, ˙ ρ ej ( v + v ′ ej ) should on average be purely radial because thephotoionized wind will be blown out preferentially along the pressure gradient. Neither F ∗ nor the viscous stresstensor is included in the momentum equation used in Paper I. A detailed comparison of the two approaches leads us toconclude that the formulation used here more properly follows the transport of energy through the turbulent cascade.After taking the second time derivative of the cloud moment of inertia, I cl = R V vir ρr dV , and substituting themomentum equation (A5) into the resulting expression, we find12 ¨ I cl = 12 ddt Z V vir ˙ ρr dV − ddt Z S vir ρ v r · d S + Z V vir r · [ ρ g + F ∗ − ∇ · (Π − T M )] dV + Z V vir ˙ ρ ej r · ( v + v ′ ej ) dV + Z V vir ˙ ρ acc r · v res dV. (A6)Upon evaluating the integrals in equation (A6) term by term, we obtain the Eulerian Virial Theorem,12 ¨ I cl = 2( T − T ) + B + W − ddt Z S vir ( ρ v r ) · d S + a I ˙ M cl R cl ˙ R cl + 12 a I ¨ M cl R + a I ˙ M ej R cl ˙ R cl + 3 − k ρ − k ρ R cl ( ˙ M ej v ′ ej − ξ ˙ M acc v esc ) (A7)Here T , T , B , and W are respectively the standard kinetic, surface kinetic, magnetic, and gravitational terms (seePaper I for precise definitions) and a I = (3 − k ρ ) / (5 − k ρ ). The final term in equation (A7) does not appear in theEVT derived in Paper I and is due to the presence of the accretion flow. This has the same form as the wind recoilterm except for the presence of a dimensionless factor ξ = Z (4 − k ρ ) x − k ρ (cid:26) f (cid:20) − x − k ρ − k ρ + Z x y ( x ′ ) x ′ dx ′ (cid:21)(cid:27) dx (A8)which arises because material is accreted at a velocity that depends on the depth of the cloud potential well. Thedimensionless variables x, y, and f are defined explicitly in Appendix B.The magnetic term B retains the form used in Paper I and derived by (McKee & Zweibel 1992) because we assumeany deformation of the magnetic field in the ambient medium caused by the presence of the cloud is negligible at thevirial surface, allowing us to approximate that the virial volume is threaded by a constant magnetic field B . Here B is the RMS value of the true magnetic field at the virial surface, which may fluctuate around B . Since material in thereservoir should also carry currents and thus generate magnetic fields, this parameterization underestimates the totalmagnetic energy. However, we expect the mean density of reservoir material within the virial volume to be an orderof magnitude smaller than the density of cloud material (see Appendix B), so the contribution of the reservoir to themagnetic energy is small compared to the contribution due to the cloud. B. PROPERTIES OF THE RESERVOIR
Consider an accreting cloud that is not forming stars and thus not generating a wind. Let M res ( r, t ) be the mass ofmaterial in the accretion flow contained within a radius r at time t and let ∆ r = v res , sys ∆ t be the distance that theaccreting gas falls in a time ∆ t . In the same time, a fraction of the material in the accretion flow will be convertedinto cloud material. In the frame comoving with the accretion flow, the change in the mass of reservoir interior to aradius r in a time ∆ t is, M res ( r, t ) − M res ( r − ∆ r, t + ∆ t ) = − ∆ t Z r ˙ ρ acc dV. (B1)Upon Taylor expanding in ∆ t and ∆ r , dropping the nonlinear terms, and evaluating the integral on the right handside of equation (B3), we find, ∆ r ∂M res ( r, t ) ∂r − ∆ t ∂M res ( r, t ) ∂t = − ˙ M acc ( r, t )∆ t. (B2)If the accretion flow is in quasi-steady state, the time derivative vanishes and M res ( r, t ) = M res ( r ). Integrating toobtain M res ( r ), we findlobal Evolution of GMCs: Accretion 21 Figure 9. (a): y (dotted line) and z (dashed line) as a function of x = r/R cl . (b): y (1) = M res ( R cl ) /M cl as a function of ζ = ˙ M acc t ff /M cl .We see M res ( R cl ) ≥ M cl for ζ & M res ( r ) = − Z r ˙ M acc ( r ′ ) v res , sys ( r ′ ) dr ′ (B3)If we wish to evaluate the above integral and obtain an expression for M res ( r ), we need to know v res , sys ( r ). Expandingequation (4), we see12 v , sys ( r ) = − Z R cl ∞ G ( M cl + M res ) r ′ dr ′ − Z rR cl GM cl ( r ′ ) r ′ "(cid:18) r ′ R cl (cid:19) − k ρ + M res ( r ′ ) M cl ( r ′ ) dr ′ . (B4)Equations (B3) and (B4) constitute a system of integral equations that must be simultaneously solved to obtain asolution for the velocity and mass profile of reservoir material within the cloud volume. If we define the functions x = r/R cl , y ( x ) = M res ( R cl x ) /M cl , and f = M cl / ( M res + M cl ), equation (B4) reduces to v res , sys ( r ) = − v esc (cid:20) f (cid:18) − x − k ρ − k ρ − Z x y ( x ′ ) x ′ dx ′ (cid:19)(cid:21) / . (B5)Defining ζ = ˙ M acc t ff /M cl and z ( x ) = R x y ( x ′ ) /x ′ dx ′ , we see that the problem of determining both v res , sys ( r ) and M res ( r ) reduces to solving the following system of nonlinear ordinary differential equations: dydx = 2 π ζx − k ρ (cid:20) f (cid:18) − x − k ρ − k ρ + z ( x ) (cid:19)(cid:21) − / , (B6) dzdx = − y ( x ) x . (B7)These equations can be solved numerically using the shooting method as follows. Since there should be no unmixedaccreted material at r = 0, we expect y (0) = 0. Since z ( x ) is defined in terms of an integral, we have no a priori knowledge of z (0). However, we do know that z (1) = 0 and we expect y (1) ≤ ζ . If we impose y (1) = c ζ where c is aconstant of order unity, we can integrate the above system and obtain a trial solution for y ( x ) and z ( x ). The constant c can then be varied until a solution is obtained with y (0) = 0. There is a unique solution for each value of ζ and thusa new solution must be obtained if ζ varies.A key assumption in this analysis was that the accretion rate is approximately constant over the cloud free-falltimescale. This is equivalent to the assumption that ζ .
1. This is a reasonable assumption, which we can see bymaking an analogy to the case of a protostellar core. Since we expect ˙ M cl ≈ M cl /t ff , r where t ff , r is the free-fall time forall of the gas in the reservoir (Stahler et al. 1980; McKee & Tan 2003) and since the reservoir should be significantly2 Goldbaum et al.more extended than the cloud, we expect the mean density of the reservoir to be much lower than the mean densityin the cloud. This implies t ff , r ≫ t ff , cl . If this condition does hold, then the condition ζ ≤ y ( x ) and z ( x ) obtained for ζ = 1 .
0. We see for this case that M res ( R cl ) ≈ . M cl , showing that even for substantial accretion rates, the gas within V cl is primarily composed ofcloud material. In Figure 9b we show y (1) = M res ( R cl ) /M cl as a function of ζ . Reservoir material becomes the primarycomponent of the volume occupied by the cloud for ζ & .
0. For this reason, we reject as unphysical any portion ofcloud history with ζ ≥
5. In practice, we find ζ . . C. DERIVATION OF THE EQUATION OF ENERGY CONSERVATION FOR AN ACCRETING CLOUD
Here we derive the global equation of energy conservation for a cloud undergoing accretion and star formation. Webegin by writing the equation of momentum conservation (equation [A5]) in Lagrangian form, ρ d v dt = − ∇ P − ρ ∇ φ + ∇ · π + J × B c + ˙ ρ ej v ′ ej + F ∗ + ˙ ρ acc ( v res − v ) . (C1)Upon contracting the above equation with v , we obtain the nonthermal energy evolution equation, ∂∂t (cid:18) ρv + ρφ cl + B − B π (cid:19) + ∇ · ρ v (cid:18) v + φ cl (cid:19) + ∇ · S P = − v · ∇ P + ρ ∂φ cl ∂t + ρ v · g res + ∇ · ( π · v ) − π : ∇ v + ˙ ρ (cid:18) v + φ cl (cid:19) + v · F ∗ + ˙ ρ ej v · v ′ ej + ˙ ρ acc ( v · v res − v ) (C2)where S P is the Poynting vector. Combining the first law of thermodynamics with the continuity equation yields theevolution equation for the thermal energy of the cloud, ∂∂t ( ρe ) + ∇ · ρ v (cid:18) e + Pρ (cid:19) = ˙ ρe + v · ∇ P + Γ − Λ (C3)where e is the internal energy per unit mass and Γ and Λ are respectively the rates of energy gain and loss per unitvolume. Here we’ve assumed that accreted material has the same thermal energy density as cloud material, which willbe true if the radiative timescales are short compared to the mechanical timescales, as in molecular cloud conditions.Summing equation (C2) and (C3), we obtain the evolution equation for the total energy of the system ∂∂t (cid:20) ρ (cid:18) v + e + 12 φ cl (cid:19) + B − B π (cid:21) + ∇ · ρ v (cid:18) v + e + Pρ + φ cl (cid:19) + ∇ · S P = 12 (cid:18) ρ ∂φ cl ∂t − ∂ρ∂t φ cl (cid:19) + ˙ ρ (cid:18) v + e + 12 φ cl (cid:19) + 12 ˙ ρφ cl + ρ v · g res + ˙ ρ ej v · v ′ ej + ˙ ρ acc ( v · v res − v ) + ∇ · ( π · v ) + v · F ∗ − Λ . (C4)Here we’ve assumed that Γ = π : ∇ v , the scalar rate of viscous dissipation. This is equivalent to the statement thatturbulent kinetic energy is converted into heat at viscous scales whereupon the energy is quickly radiated away. Whilethere may be other heating mechanisms, including cosmic rays and protostellar radiation, we neglect these sources bynoting that any local heating should be offset in a time much shorter than the dynamical time.If we define the total energy due to the presence of the cloud, E cl = Z V cl ρ (cid:18) v + e + 12 φ cl (cid:19) dV + Z V vir B − B π dV (C5)and the total energy in the ambient medium, E amb = Z V vir − V cl ρ (cid:18) v + e + 12 φ cl (cid:19) dV, (C6)we can integrate equation (C4) over the virial volume and obtain an evolution equation for E cl . Portions of thiscalculation are performed explicitly in Paper I and we do not reproduce those results here. The new terms stemfrom our inclusion of accretion as well as the inclusion of F ∗ , the random component of the stellar forcing term. Theintegrals over the new terms are evaluated below.The first new term is due to the gravitational influence of the reservoir, Z V vir ρ v · g res dV = − χ ˙ R cl R cl GM R cl (C7)wherelobal Evolution of GMCs: Accretion 23 χ = (3 − k ρ ) Z x − k ρ y ( x ) dx. (C8)This is the gravitational work done on the cloud by the reservoir material as the cloud expands and contracts.The random stellar forcing term can be evaluated by noting that r and F ∗ are assumed to be uncorrelated, so R V vir v · F ∗ dV = R V vir v turb · F ∗ dV . This integral depends on the degree of correlation between two randomly orientedvectors, v turb and F ∗ . Since, at a fixed time, the direction of expansion of an H ii region is the same as the direction ofmomentum injection, these vectors should indeed be highly correlated. This term then corresponds to the net injectionof energy by H ii regions, G cl = Z V vir v turb · F ∗ dV. (C9)For the terms proportional to ˙ ρ acc , we have Z V vir ˙ ρ acc v · v res dV = − (cid:18) − k ρ − k ρ (cid:19) ξ ˙ R cl ˙ M acc v esc + Z V vir ˙ ρ acc v turb · v res dV (C10)and Z V vir ˙ ρ acc v dV = a I ˙ M acc ˙ R cl + 3 ˙ M acc σ (C11)where we’ve used our assumption that v res , rand is randomly oriented with respect to r .We are left with one last integral on the right hand side of equation (C10) that depends on the correlation betweentwo vector fields, v turb and v res . Suppose v res and v turb are perfectly correlated. In that case, the integral we areconcerned with can be readily evaluated in three limits, | v res | = | v turb | , | v res | ≫ | v turb | , and | v res | ≪ | v turb | . In thefirst limit | v res | = | v turb | the transfer of material to the cloud is merely an act of relabeling, so v turb · v res = v . Nowconsider the limit | v res | ≫ | v turb | . We approximate that a parcel of reservoir material mixes with the cloud once it hasswept up a mass of cloud material equal to its own mass. Since the interaction between cloud material and reservoirmaterial must be inelastic if they are to mix, the velocity of cloud material must be driven by the act of mixing to v turb = v res /
2. Conversely, if | v res | ≪ | v turb | , the reservoir material is driven to a velocity v res = v turb /
2. We obtainthe correct answer in all three limits if we assume Z V vir ˙ ρ acc v turb · v res dV = 12 Z V vir ˙ ρ acc ( v + v ) dV. (C12)This is an upper bound on the value of the integral in the limit of perfect correlation between v turb and v res . If v turb and v res are uncorrelated then the integrand should on average be zero. Thus we have,0 ≤ Z V vir ˙ ρ acc v turb · v res dV ≤ Z V vir ˙ ρ acc ( v + v ) dV. (C13)To evaluate this integral, we linearly interpolate between the upper limit and lower limit, Z V vir ˙ ρ acc v turb · v res dV = ϕ Z V vir ˙ ρ acc ( v + v ) dV = ϕ ( 32 ˙ M acc σ + 32 ˙ M acc σ + γ ˙ M acc v ) (C14)where ϕ is an interpolation parameter that ranges between zero and unity and γ = 12 + f (cid:18) − k ρ − χ − k ρ (cid:19) . (C15)Summing the individual terms derived above yields the global form of the equation of energy conservation, d E cl dt = ˙ M cl M cl [ E cl + (1 − η B ) W ] + GM cl ˙ M cl R cl χ − M cl ˙ R cl ˙ M cl R cl ! + (cid:18) − k ρ − k ρ (cid:19) ˙ R cl ( ˙ M ej v ′ ej − ξ ˙ M acc v esc ) − πP amb R ˙ R cl − a I ˙ M acc ˙ R − M acc σ cl + ϕ ( 32 ˙ M acc σ + 32 ˙ M acc σ + γ ˙ M acc v ) + G cl − L cl (C16)where L cl = R V vir Λ dVdV