ISM Structure & Growth Histories

# An Excursion-Set Model for the Structure of GMCs and the ISM

## Abstract

The ISM is governed by supersonic turbulence on a range of scales. We use this simple fact to develop a rigorous excursion-set model for the formation, structure, and time evolution of dense gas structures (e.g. GMCs, massive clumps and cores). Supersonic turbulence drives the density distribution in non-self gravitating regions to a lognormal with dispersion increasing with Mach number. We generalize this to include scales (the disk scale height), and use it to construct the statistical properties of the density field smoothed on a scale . We then compare conditions for self-gravitating collapse including thermal, turbulent, and rotational (disk shear) support (reducing to the Jeans/Toomre criterion on small/large scales). We show that this becomes a well-defined barrier crossing problem. As such, an exact “bound object mass function” can be derived, from scales of the sonic length to well above the disk Jeans mass. This agrees remarkably well with observed GMC mass functions in the MW and other galaxies, with the only inputs being the total mass and size of the galaxies (to normalize the model). This explains the cutoff of the mass function and its power-law slope (close to, but slightly shallower than, ). The model also predicts the linewidth-size and size-mass relations of clouds and the dependence of residuals from these relations on mean surface density/pressure, in excellent agreement with observations. We use this to predict the spatial correlation function/clustering of clouds and, by extension, star clusters; these also agree well with observations. We predict the size/mass function of “bubbles” or “holes” in the ISM, and show this can account for the observed HI hole distribution without requiring any local feedback/heating sources. We generalize the model to construct time-dependent “merger/fragmentation trees” which can be used to follow cloud evolution and construct semi-analytic models for the ISM, GMCs, and star-forming populations. We provide explicit recipes to construct these trees. We use a simple example to show that, if clouds are not destroyed in crossing times, then all the ISM mass would be trapped in collapsing objects even if the large-scale turbulent cascade were maintained.

###### keywords:
galaxies: formation — star formation: general — galaxies: evolution — galaxies: active — cosmology: theory

## 1 Introduction

The origins and nature of structure in the interstellar medium (ISM) and giant molecular clouds (GMCs) represents one of the most important unresolved topics in both the study of star formation and galaxy formation. In recent years, there have been several major advances in our understanding of the relevant processes. It is clear that a large fraction of the mass in the ISM is supersonically turbulent over a wide range of scales, from the sonic length (pc) through and above the disk scale height (kpc). A generic consequence of this super-sonic turbulence – so long as it can be maintained – is that the density distribution converges towards a lognormal PDF, with a dispersion that scales weakly with mach number (e.g. Vazquez-Semadeni, 1994; Padoan et al., 1997; Scalo et al., 1998; Ostriker et al., 1999).

Without continuous energy injection, this turbulence would dissipate in a single crossing time, and the processes that “pump” turbulence (generally assumed to be related to feedback from massive stars) remain poorly understood (see e.g. Mac Low & Klessen, 2004; McKee & Ostriker, 2007; Hopkins et al., 2012, and references therein). However, provided this turbulence can be maintained, it is able to explain the relatively small fraction of mass which collapses under the runaway effects of self-gravity and cooling (Vázquez-Semadeni et al., 2003; Li et al., 2004; Li & Nakamura, 2006; Padoan & Nordlund, 2011). In this picture, star formation occurs within dense cores, themselves typically embedded inside giant molecular clouds (GMCs), which represent regions where turbulent density fluctuations become sufficiently overdense so as to be marginally self-gravitating and collapse (Evans, 1999; Gao & Solomon, 2004; Bussmann et al., 2008). Some other process such as stellar feedback is believed to be responsible for disrupting the clouds, after a few crossing times (e.g. Evans et al., 2009). The turbulent cascade has also been invoked to explain GMC scaling relations, such as the size-mass and linewidth-size relations (Larson, 1981; Scoville et al., 1987).

However, despite this progress, there remains no rigorous analytic theory that can simultaneously predict these properties, as well as other key observables such as the GMC mass function, and the spatial distribution of gas over and under-densities in the ISM.

The approximately Gaussian distribution of the logarithmic density field, though, suggests that considerable progress might be made by adapting the excursion set or “extended Press-Schechter” formalism. This has proven to be an extremely powerful tool in the study of cosmology and galaxy evolution. The seminal work by Press & Schechter (1974) derived the form of the halo mass function via a simple (albeit somewhat ad hoc) calculation of the mass fraction expected to be above a given threshold for collapse, expected in a Gaussian overdensity distribution with the variance as a function of scale derived from the density power spectrum. Bond et al. (1991) developed a rigorous analytic (and statistical Monte Carlo) formulation of this, defining the excursion set formalism for dark matter halos. Famously, this resolved the “cloud in cloud” problem, providing a means to calculate whether structures were embedded in larger collapsing regions. Since then, excursion set models of dark matter have been studied extensively: they have been generalized and used to predict – in addition to the halo mass function – the spatial distribution/correlation function of halos (Mo & White, 1996), the distribution of voids (Sheth & van de Weygaert, 2004), the evolution and structure of HII regions in re-ionization (Haiman et al., 2000; Furlanetto et al., 2004), and many higher-order properties used as cosmological probes. By incorporating the time-dependence of the field, they have been used to study the growth and merger histories of halos and to construct Monte Carlo “merger trees” (Bower, 1991; Lacey & Cole, 1993). These trees formed the basis for the extensive field of semi-analytic models for galaxy formation, in which analytic physical prescriptions for galaxy evolution are “painted onto” the background halo evolution (e.g. Somerville & Kolatt, 1999; Cole et al., 2000). It is not an exaggeration to say that it has proven to be one of the most powerful theoretical tools in the study of large scale structure and galaxy formation.

There have been other, growing suggestions of similarities between the mathematical structure of the ISM and that invoked in excursion set theory. The mass function of GMCs, for example, has a faint-end slope quite similar to that of galaxy halos (both close to ), suggestive of hierarchical collapse. Vazquez-Semadeni (1994) attempted to rigorously examine whether the structure of the ISM should be “hierarchical,” although they strictly define this as the probability of many independent fluctuations dominating the “peaks” in the density distribution (which does not technically need to be satisfied in excursion set theory). This is related to (but not equivalent to) the large body of work on the quasi-fractal structure of the ISM (see e.g. Elmegreen, 2002, and references therein). On smaller scales, Krumholz & McKee (2005) suggested that the fraction of a lognormal PDF above a “collapse threshold” at the sonic length could explain the fractional mass forming stars per free-fall time, inside of GMCs. Padoan et al. (1997); Padoan & Nordlund (2002) suggested that the distribution of lognormal density fluctuations above a threshold overdensity could explain the shape of the stellar initial mass function (IMF). Scalo et al. (1998) explicitly discuss the analogy between this and cosmological density fluctuations, and Hennebelle & Chabrier (2008) expanded upon the Padoan et al. (1997) argument using a derivation almost exactly analogous to the original Press & Schechter (1974) derivation, and showed that it agreed well with the standard IMF.

But despite these suggestions, and the enormous successes of the excursion set model in cosmological applications, there has been no attempt to translate the excursion set formalism to the problem of the ISM and GMC evolution. At first glance, it is obvious why. The cosmological excursion set theory is applied to small fluctuations of the linear density field, in the linear regime, to dark matter (collisionless) systems, with Gaussian, nearly scale-free fluctuations seeded by inflation, and to Lagrangian “halos” which (modulo mergers) are conserved in time. The Gaussian distribution of ISM densities represent large fluctuations in the logarithmic density field, which are a product of a fully non-linear, turbulent, gaseous (collisional) medium, and evolve both rapidly and stochastically in time.

However, in this paper, we will show that although the physics involved are very different, none of these differences fundamentally invalidates the underlying mathematical formalism of the excursion set theory.

Here, we develop a rigorous excursion-set model for the formation, structure, and time evolution of structures in the ISM and within GMCs. We show that this is possible, and that it allows us to develop statistical predictions of ISM properties in a manner analogous to the predictions for the halo mass function. In § 2 we describe the model. First (§ 2.1), we derive the conditions for self-gravitating collapse in a turbulent medium (the “collapse threshold”), in a manner generalized to both small (sonic length) and large (above the disk scale height) scales. Next (§ 2.2), we discuss the density field, and, assuming it has a lognormal character, construct the statistical properties of the field smoothed on a physical scale , which allows us to define the excursion set “barrier crossing” problem. In § 3, we use this to derive an exact “self-gravitating object” mass function, over the entire range of masses (from the sonic length to disk mass), and show that it agrees remarkably well with observed GMC mass functions and depends only very weakly on the exact turbulent properties of the medium (including deviations from a lognormal PDF). In § 4, we show that the model also predicts the linewidth-size and size-mass relations of GMCs, and their dependence on external galaxy properties. We also examine how this depends on the exact properties of the turbulent cascade. In § 5, we extend the model to predict the spatial correlation function and clustering properties of clouds (and, by extension, young star clusters), and compare this to observations. In § 6, we predict the size and mass distributions of underdense “bubbles” or “holes” in the ISM which result simply from the same normal turbulent motions. We show that this can explain most or all of the distribution of HI “holes” observed in nearby galaxies, without explicitly requiring any feedback mechanism to power the hole expansion. In § 7, we generalize the model to construct time-dependent “GMC merger/fragmentation trees” which follow the time evolution, growth histories, fragmentation, and mergers of clouds. In § 7.2 we provide simple recipes to construct these trees, and discuss how they can be used to build semi-analytic models for GMC and ISM evolution and star formation, in direct analogy to semi-analytic models for galaxy formation. We use a very simple example of this to predict the rate at which the gas in the ISM collapses (absent feedback) into bound structures, show that this agrees well with the results of fully non-linear turbulent box simulations, and argue that feedback must destroy clouds on a short timescale (a few crossing times) to prevent runaway gas consumption. Finally, in § 8, we summarize our results and conclusions and discuss a number of possibilities for future work, both to improve the accuracy of these models and to enable predictions for additional properties of the ISM.

## 2 The Model

The fundamental assumption of our model is that non-rotational velocities are dominated by super-sonic turbulence (down to some sonic length), with some power spectrum or 1 which is maintained by any process (presumably stellar feedback) in approximate statistical steady-state. As we discuss in § 8, all other assumptions we make are convenient approximations to simplify our calculations, but it is possible to generalize the model.

The two key quantities we need to calculate the cloud mass function and other properties are the conditions for “collapse” of a cloud (i.e. conditions under which self-gravity can overcome turbulent forcing) and the power spectrum of density fluctuations. Below, we show how these can be calculated for a turbulent medium from the velocity power spectrum; however, in principle they can be completely arbitrary (for example, specified ad hoc from numerical simulations or observations). So long as they are known, the rest of our model proceeds identically.

### 2.1 Collapse in a Turbulent Medium

First, for simplicity, consider gas in a galaxy whose average properties are evaluated on a scale where the velocity dispersion is highly supersonic (, where is the sonic length), but where shear from the disk rotation and large-scale density gradients can be neglected (, where is the disk scale-height). The turbulent dispersion on these scales is . If the turbulence has a power-law cascade over this interval then . If the region has some mean density (on the same scale ), then the potential from self-gravity is while the kinetic energy in turbulence is ; the region will be gravitationally bound and “trapped” when , where depends on the shape (internal structure) of the density perturbation. Formally, we also need to check whether the momentum “input” rate from the turbulent cascade (equal to the dissipation rate in steady-state) is less than the gravitational force, and whether the energy input rate is less than the rate at which a gravitationally collapsing object will dissipate. However, because for super-sonic turbulence, the timescale for energy or momentum dissipation on a scale just scales with the crossing time , we obtain the identical dimensional scaling for all of these criteria.2

These are simply a restatement of the Jeans criterion, for wavenumber , but with the sound speed replaced by the turbulent velocity dispersion . For an individual -mode (sinusoidal density perturbation), the criteria becomes

 ρ(R)≥k2⟨v2t(k)⟩4πG∝k3−p∝Rp−3 (1)

where the latter equalities assume a power-law spectrum (Vazquez-Semadeni & Gazol, 1995). If the system is marginally stable with density on scale , then this simply becomes . If we are in the super-sonic regime, then we expect something like Burgers turbulence (Burgers, 1973), with ; but we will discuss this further below.

Now generalize this to a more broad range of radii. On small scales, we need to include the effects of thermal pressure: this amounts to a straightforward modification of the Jeans criterion with (Chandrasekhar, 1951; Bonazzola et al., 1987).3 On large scales, we need to include the effects of rotation stabilizing perturbations. If we focused only on very large () scales, where we can neglect the disk thickness, then we simply re-derive the (Toomre, 1977) dispersion relation and collapse conditions, with the gas “dispersion” . More generally, Begelman & Shlosman (2009) note that the dispersion relation for growth of density perturbations in a turbulent disk (with finite thickness ) can be written:

 ω2 =κ2+σg(k)2k2−2πGΣ|k|1+|k|h (2) =κ2+σg(k)2k2−4πGρ|k|h1+|k|h (3)

where is the disk surface density, the average density on scale , and the disk scale-height, the turbulent velocity dispersion, and the usual epicyclic frequency. This differs from the infinitely thin-disk dispersion relation by the term , which accounts for the finite scale height for modes with scales (Vandervoort, 1970; Elmegreen, 1987; Romeo, 1992).4 Note that this relation nicely interpolates between the Jeans criterion which we derived above on small scales (), and the Toomre (thin-disk) dispersion relation on large scales ().

If the average density is , and corresponding average surface density , then we can define the usual Toomre at the scale

 Q0(h)≡κσg(h)πGΣ0=h~κΩ2πGΣ0 (4)

where the second equality follows from , which is be true for any disk in vertical equilibrium, and we define ( for a constant- disk). If we define the convenient dimensionless form of , , we can write the criterion for instability () as

 ρρ0 ≥ρcρ0≡Q0(h)2~κ(1+~k)[σg(k)2σg(h)2~k+~κ2~k−1] (5)

Note that the assumption of a finite ensures that so long as there is any non-gaseous component of the potential, the gas alone is not self-gravitating on arbitrarily large scales (this is important below, to un-ambiguously define the largest self-gravitating scales of clouds). Again, on small scales , this reduces to the Jeans criterion , and on large scales it becomes .

Kim et al. (2002) note that is straightforward to further generalize this criterion to include the effects of magnetic fields by taking , where is the Alfvèn speed. If we follow the usual convention in the literature and assume is constant, then changing the strength of magnetic fields is identically equivalent to changing the sound speed/mach number (which we explicitly consider below). Even if we allow to have an arbitrary power spectrum, the results are quite similar to this renormalization – for any power spectrum where the magnetic energy density is peaked on large scales, it is nearly equivalent to renormalizing the turbulent velocities; for a power spectrum peaked on small scales, equivalent to renormalizing the sound speed. We therefore will not explicitly consider magnetic fields in what follows, but emphasize that they are straightforward to include if their power spectrum is known.

Formally, the turbulent velocity power spectrum must eventually flatten/turn over on large scales , both by definition (since itself traces the maximal three-dimensional dispersions) and to avoid energy divergences. If it did not, we would recover on large scales in gas-rich systems! Constancy of energy transfer and energy conservation require that the slope become at least as shallow as . A good approximation to the behavior seen in simulations is obtained by generalizing the exact correction for near the lowest wavenumbers in the inertial scale in Kolmogorov turbulence (Bowman, 1996), taking , which interpolates between these regimes. This may not be exact. Fortunately however, even if we ignored this correction entirely, we can see immediately from Equation 5 that for any reasonable power spectrum (), the dominant velocity/pressure term on scales is the disk shear (), not . We therefore include this turnover, but stress that it is not necessary to our derivation and has only weak effects on our conclusions.

### 2.2 The Density Distribution

The other required ingredient for our model is an estimate of the density PDF/power spectrum. We emphasize that the our methodology is robust to the choice of an arbitrary PDF and/or power spectrum in . We could, for example, simply extract a density power spectrum (or fit to it) from simulations or observations. This is, however, less predictive – so in this paper, we chose to focus on the case of supersonic turbulence in which case it is possible to (at least approximately) construct the density PDF knowing only the velocity power spectrum information.

As discussed in § 1, in idealized simulations of supersonic turbulence with a well-defined mean density and mach number on a scale , the distribution of densities tends towards a lognormal distribution

 dp(δ|k) =1σk√2πexp(−δ22σ2k)dδ (6) δ (7)

where because is the mean density,

 (8)

This form of the PDF and our results are identical whether we define all quantities as volume-weighted or mass-weighted, so long as we are consistent throughout: here it is convenient to define all properties as volume-weighted (otherwise is scale-dependent).

The dispersion in these simulations is a function of the rms (one-dimensional) Mach number averaged on the same scale ,

 σk≈(ln[1+34M(k)2])1/2 (9)

which is naturally expected for supersonic turbulence with efficient cooling (because the variance in in “events” – namely strong shocks – scales as ).5

If the turbulence obeys locality – i.e. if the density distribution averaged on some small scale depends only on the local gas properties on that scale as opposed to e.g. the structure on much larger scales – then the distribution of densities averaged over any spatial scale with some window function is also a lognormal in , with variance

 σ2(R)=∫dln(k)σ2k(M[k])|W(k,R)|2 (10)

where is the Fourier transform of . This is easy to see if we recursively divide an initially large volume (e.g. the entire disk) into sub-regions with different mean and on scale ; each of these sub-regions is a “box” that should obey the density distributions above, and so on. Because it greatly simplifies the algebra, we will generally follow the standard practice in the excursion set literature and choose to be a Fourier-space top-hat: if and if . This choice is arbitrary, but so long as it is treated consistently, our subsequent results are essentially identical (we will show, for example, that using a Gaussian window function makes a small difference in all predicted quantities).6

It should immediately be clear, however, that if we simply extrapolated , the dispersion would be divergent! Physically, this would imply ever larger fluctuations in on arbitrarily large scales; but this cannot be true once the scale approaches that of the entire disk. As , the fact that the disk has finite mass means that . The resolution of this apparent dilemma is evident in Equation 2: what matters in in the dispersion is the effective “pressure” from ; on sufficiently large scales , the differential rotation plays an identical physical role. We can therefore generalize Equation 9 as

 σk ≈(ln[1+34v2t(k)c2s+κ2k−2])1/2 (11) =(ln[1+34M2(k)1+M2h~κ2/|kh|2])1/2 (12)

where . This ensures the correct physical behavior, as for all plausible turbulent .

At some level, our assumptions must break down. And although it is well-established that the density PDF at the resolution limit in numerical simulations (in a “box-averaged” sense) approaches the behavior of Eqns. 6-9, it is less clear whether we can assume this on a -by- basis and so derive Eqns. 10-11. The lognormal character of the density distribution holding on various smoothing scales as we assume is, however, supported in the investigations of Lemaster & Stone (2009) and Passot & Vazquez-Semadeni (1998); Scalo et al. (1998). And any distribution which is lognormal in either real space or space must be lognormal in both. Moreover, the robustness of this assumption is supported by the conservation of lognormality in resolution studies, since all simulations essentially measure the PDF smoothed over a window function corresponding to their resolution limits. To the extent that there is some violation of these assumptions in e.g. the higher-order-structure functions (although they are largely consistent with locality when is large; see Boldyrev et al., 2002; Padoan et al., 2004; Schmidt et al., 2008), this is really a question of the degree to which the density PDF globally departs from a lognormal, which we discuss below.

What is somewhat less clear is how Eq. 9 generalizes on a scale-by-scale basis. Analytically, the same arguments that prove the density distribution of isothermal turbulence should converge to a lognormal with real-space variance trivially generalize to a -by- basis (Eq. 9; see Passot & Vazquez-Semadeni, 1998; Nordlund & Padoan, 1999). If locality also holds, Eq. 10 must follow. This is the origin of the expectation for analytic models of the density power spectrum. Note that, as defined, is equivalent to the logarithmic density power spectrum, . When is not large, in Eq. 9 scales , so . This is just the well-known expectation that in the weakly compressible regime, the log density power spectrum should have the same shape as the velocity power spectrum. Kowal et al. (2007) and Schmidt et al. (2009) show that this is a good approximation for the field in simulations of supersonic turbulence. At higher , this should flatten logarithmically, and this is seen in numerical simulations in Kowal et al. (2007), in excellent agreement with Eq. 9. These behaviors, and the approximate normality of , appear to hold even in simulations which include explicitly non-local effects such as magnetic fields, self-gravity (excluding the collapsing regions), radiation pressure, photoionization, and non-isothermal gas with realistic heating/cooling (see e.g. Ostriker et al., 1999; Klessen, 2000; Lemaster & Stone, 2009; Hopkins et al., 2012).

Even if our analytic derivation is not exact, we can think of the resulting and implied log density power spectrum () as a convenient approximation for the power spectrum measured in hydrodynamic simulations and observations. At sufficiently large , where is small, ; a steep falloff with for typical ; at smaller (but still smaller scales than the disk scale-height) is large and this flattens to with a small logarithmic correction. This is exactly the behavior directly measured in numerical simulations (Kowal et al., 2007; Schmidt et al., 2009). Qualitatively similar behavior is seen in the linear density spectrum, but it is important to distinguish the two, since it is well known that large fluctuations at higher will further flatten the linear spectrum (see Scalo et al., 1998; Vázquez-Semadeni & García, 2001; Kim & Ryu, 2005; Kritsuk et al., 2007; Bournaud et al., 2010). It is also consistent with observations of the projected surface density power spectrum in local galaxies and star-forming regions (Stanimirovic et al., 1999; Padoan et al., 2006; Block et al., 2010). If we integrate to get we obtain constant as , with an absolute value of  dex for a range and . This range is quite similar to the range measured in on the smallest resolved scales in a wide range of simulations that have a sufficiently large dynamic range in scales to probe the typical Mach numbers in GMCs and disk scale heights (see Vazquez-Semadeni, 1994; Nordlund & Padoan, 1999; Ostriker et al., 2001; Mac Low & Klessen, 2004; Slyz et al., 2005; Hopkins et al., 2012). It also agrees well with measured values of the dispersion in the real ISM (Wong et al., 2008; Goodman et al., 2009a; Federrath et al., 2010).

## 3 The Mass Function

The question of the mass collapsed on different scales is now a well-posed barrier crossing problem. The quantity – the logarithm of the density smoothed on the scale – is distributed as a Gaussian random field with variance and zero mean, with a well-defined barrier

 δc(R)≡δ(ρc,R)=ln(ρcρ0)−⟨ln(ρρ0)⟩ (13)

which, upon crossing, leads to collapse. The mass of a collapsed object is simply the integral of the density over the effective volume of a window of effective radius in real space. If the medium were infinite and homogenous, this would just be ; however, we need to account for the finite vertical thickness of the disk. For the same vertical exponential profile that gives rise to the dispersion relation in Equation 2, the total mass inside is

 M(ρ|Rw)≡4πρ(Rw)h3[R2w2h2+(1+Rwh)exp(−Rwh)−1] (14)

where is the midplane density (chosen for consistency with the dispersion relation). This formula simply interpolates between for and for , as it should.

The fraction of the total mass which is in collapsed objects, averaged over a given smoothing scale , is then just

 Fcoll(Rw)=1Mtot∫∞δcM(ρ|Rw)p(δ|Rw)dδ (15)

where . Naively, we would equate this to the mass function of such objects with the relation . Indeed, up to a normalization factor, that is exactly the original approach of Press & Schechter (1974). However, this neglects the “cloud in cloud” problem: namely, it does not resolve whether or not a collapsed region on a scale is independent, or is simply a random sub-region of a larger object collapsed on a scale . For the case of a constant , accounting for this amounts to a simple re-normalization; but there is no simple closed-form analytic solution for the complicated here, and we will show that accounting for this behavior is critical.

### 3.1 Exact Solution

To derive the exact mass function solution we turn to the standard Monte Carlo excursion-set approach. Consider the density field at some arbitrary location , smoothed over some window corresponding to the radius (and mass ) . This is the convolution ; so if we Fourier transform, we obtain . In other words, the amplitude is simply the integral of the contribution from all Fourier modes , weighted by the Fourier-space window function.

In this sense, we can think of the (statistical) evaluation of the density field as the results of a “random walk” through Fourier space. Bond et al. (1991) show that this integration becomes particularly simple for the case of a Gaussian field with a Fourier-space tophat window, in which case the probability of a transition from to as we step from a scale to is given by

 p(δ1+Δδ)dΔδ =1√2πΔSexp(−(Δδ)22ΔS)d(Δδ) (16) ΔS ≡S2−S1≡σ2(R2)−σ2(R1) (17)

where we define the variance

 S(R)≡σ2(R) (18)

i.e. the increment is a Gaussian random variable with standard deviation .

If we begin on some sufficiently large initial scale (), then the overdensity and density variations must go to zero. We then have the well-defined initial conditions for the walk, , . Starting at some arbitrarily large , and moving to progressively smaller scales with increments7 in or ( or ) we can then compute the trajectory or ,

 δ(Ri)≡Rj>Ri∑jΔδj . (19)

At each scale , we then evaluate whether or not the barrier has been crossed,

 δ(Ri)≥δc(Ri) . (20)

If this is satisfied, we then associate that trajectory with a collapse on the scale and mass .

Recall, we are sampling the field , so the fraction of trajectories that cross the barrier in some interval or (equivalently) represents the probability of an Eulerian volume element being collapsed on that scale. This corresponds to a differential mass . Since the total mass associated with the mass function is , we have the predicted mass function or “first-crossing distribution”:

 dndM=ρc(M)MdfdM (21)

where is the differential fraction of trajectories that cross between and .

This formalism has several advantages. It provides an exact solution that also allows us to rigorously calculate the normalization and shape of the mass function. It also allows us to explicitly resolve the “cloud-in-cloud” problem, i.e. to address the situation where a trajectory crosses the barrier multiple times. Figure 1 plots the resulting mass function (for a few choices of parameters, which just determine the normalization of the mass function and will be discussed below). We also compare the mass function if we were to ignore the “cloud in cloud” problem – i.e. where we treat every crossing above on a smoothing scale as a separate cloud. At the highest masses, the difference is small – this is because the variance is small and is large, so the probability of being inside a “yet larger” cloud vanishes. However, at lower masses, the difference rapidly becomes quite large (order-of-magnitude) – much larger than the factor of the Press-Schechter mass function. This owes to the complicated behavior of , which increases again on small scales. Failure to properly account for the cloud-in-cloud problem and moving barrier will clearly lead to large inaccuracies.

### 3.2 Key Behaviors

If the barrier were constant, the mass function of collapsed objects would then simply follow the Press-Schechter formula;

 dnPSdM=ρ0M2√2πδcσ∣∣dlnσdlnM∣∣exp(−ν22) (22)

where is the collapse threshold in units of the standard deviation () of the smoothed density field on the scale corresponding to .

However, the barrier here is not constant (it depends on ). A reasonable approximation to the first-crossing distribution, however, is given by

 dndM ≈ρcM2√2π~Bσ∣∣dlnσdlnM∣∣exp(−ν22) (23) ~B ≡{ln(ρc,min/ρ0)M

where is the critical density at the most-unstable scale. This is motivated by the exact solution for the first-crossing distribution for a linear barrier with , but with held constant below .8 Because the deviation from a constant barrier is only logarithmic, these formulae do not differ too severely, and we can gain considerable insight from their functional forms.

Consider the behavior of both and , which define three primary regimes. On scales above the sonic length but below , most of the dynamic range for GMCs, (for power-law turbulent cascades) is large, so is a very weak function of (most of the contribution comes from the largest scales, since ) while decreases with so . Therefore is a (logarithmically) decreasing function of mass. So we expect an approximately power-law mass function with slope . This implies similar mass per logarithmic interval in mass and simply follows from gravity – which is self-similar – being the dominant force (since the turbulence is super-sonic). To the extent that the slope deviates from , it is because the barrier gets larger towards lower . From the above equation, with

 α ≈−2+ln(M0/M)(3−p)2/2p2σ2 (28) ≈−2+0.1log(M0/M) (29)

(where is approximately the location of the mass function “break”; formally for MW-like systems). In other words, we expect a slope which is shallower than by a small logarithmic correction, , as observed.

At very small scales we approach the sonic length, ; the growth in becomes vanishingly small () while continues to increase logarithmically as before. The mass function must therefore flatten or turn over, with a rapidly decreasing mass in clouds below the sonic length (although the absolute number may still rise weakly).

At large scales above , decreases rapidly with increasing – the contribution from large scales goes as as , while now also increases (so ), so the mass function is exponentially cut off as . We caution that at the largest size/mass scales, global gradients in galaxy properties – which are currently neglected in our derivation of the collapse criterion – may become significant. However, the number of clouds in this limit is small.

### 3.3 Comparison with Observations

Figure 1 plots the predicted mass function: we show the exact solution, both excluding and including “clouds in clouds,” and the approximations in Equation 2327. For our “standard” model, we will assume the disk is marginally stable (), and that the turbulence, being supersonic and rapidly cooling, should have (see the discussion in § 1). Motivated by observations, we normalize the turbulent spectrum by assuming a Mach number on large scales (though we will show this exact choice has very weak effects, provided ). With these choices, the model is completely fixed in dimensionless terms. To predict an absolute number and mass scale of the mass function, we require some normalization for the galaxy properties: some measure of the local gas properties (mean density, velocity dispersion, surface density, etc, to set the mass and spatial scales) and total galaxy mass or size (to know the gas mass available). Because of our assumption of marginal stability, many of these properties are implicitly related – we need only specify e.g. a total disk mass, gas fraction, and spatial size. Or equivalently, a mean density, velocity dispersion, and total mass.

Taking typical observed values for the total gas mass, mean density, and velocity dispersion in the Milky Way, we plot the resulting predicted GMC mass function and compare to that observed. Because we are considering the total gas mass of the inner MW, we need to compare with a GMC mass function corrected to the same effective volume – we therefore compare with the values in Williams & McKee (1997) (who attempt to construct a “galaxy-wide” GMC mass function for the same total volume). We then repeat the experiment with the average properties observed in the LMC and M33, and compare with the mass function compilations in Rosolowsky (2005); Fukui et al. (2008), corrected to the appropriate survey area.

In each case, the predicted mass functions agree remarkably well with the observations. We emphasize that although the observed densities and masses enter into the normalization of the mass function, the shape, which agrees extremely well, is entirely an a priori prediction. Moreover, the assumed densities do not entirely determine the normalization – because the barrier and variance are finite at all radii, the models here specifically predict that not all mass is in bound units. In fact, only of the total mass is predicted to be in such units – for the MW, the total bound GMC mass is predicted to be , in good agreement with that observed (Williams & McKee, 1997). Likewise, the details of our stability and collapse conditions determine where, relative to the Jeans mass, the “break” in the mass function occurs.9

We should caution that it is not entirely obvious that our predicted mass function is the same as that observed. The mass function here is well-defined because we restrict to self-gravitating objects and resolve the cloud-in-cloud problem, knowing the three-dimensional field behavior (and assuming spherical collapse). In the observations, the methods used to distinguish substructures and the choice of how to average densities (in spherical or arbitrarily shaped apertures) can make non-trivial differences to the mass function (Pineda et al., 2009). This may be considerably improved by the use of more sophisticated observational techniques that attempt to statistically identify only self-gravitating structures (see e.g. Rosolowsky et al., 2008); preliminary comparison of these methods in hydrodynamic simulations and observations suggests that most of the identified GMCs are indeed self-gravitating structures so the key characteristics of the GMC MFs in our comparison should be robust, although details of individual clouds may change significantly (Goodman et al., 2009b).

### 3.4 Effects of Varying Assumptions

Of course, it is important to check how sensitive the predicted mass functions are to the assumptions in our model. Figure 2 shows the results of varying these assumptions. We plot the mass function in dimensionless units (, with the absolute mass being an arbitrary normalization).

If we assume Kolmogorov turbulence ( instead of ), the predicted mass function is nearly identical at intermediate and high masses, but flattens more rapidly at low masses, because the velocity drops more slowly at small scales so rises more steeply. The difference agrees well with the scaling in Equation 28, which predicts a faint end slope for instead of for .

If we vary the Mach number on large scales (or, equivalently, the assumed sound speed or magnetic field strength), the differences are very small at almost all masses, because the assumption that the disk as a whole is marginally stable effectively scales out the absolute value of . What does determine is the (dimensionless) scale of the sonic length (), below which the mass function will flatten. With lower , this happens at higher masses – but still quite low in absolute terms (, or  dex below the maximum GMC masses).

As noted above, the exact manner in which the velocity power spectrum should flatten at large scales is uncertain. We therefore re-calculate the mass function ignoring such flattening entirely – i.e. assuming for all . This makes the very high-mass end of the mass function slightly more shallow, but has a negligible effect at all other masses. Since the only difference will be in the regime where the number of clouds is (so subject to large Poisson fluctuations) it is difficult to constrain this from observations.

Re-calculating our results with a different window function makes little difference. We test this with a Gaussian window function (convenient as it remains a Gaussian in real and Fourier space). As discussed in Zentner (2007), this makes the calculation more complex because we can no longer treat the Fourier-space trajectory as having uncorrelated steps; following Bond et al. (1991) the first-crossing distribution is computed by numerically integrating a Langevin equation. However, we hold our mass definition fixed; with this choice, for fixed , the exact choice of window shape about introduces only small () corrections (we refer to the discussion therein and Maggiore & Riotto 2010a for more detailed discussion of the effects of different window functions).

What if the density distribution is not a lognormal? It has been suggested, for example, that for systems which have significant gas pressure and whose equations of state are non-isothermal, or which have large magnetic fields, the density distribution may more closely resemble a power-law (see e.g. Passot & Vazquez-Semadeni, 1998; Scalo et al., 1998; Ballesteros-Paredes et al., 2011b). This is certainly still treatable with the excursion set formalism: there has been considerable discussion in the literature regarding the halo mass function and bias with non-Gaussian primordial fluctuations (see Matarrese et al., 2000; Afshordi & Tolley, 2008; Maggiore & Riotto, 2010b, and references therein). However most of these rigorous approaches assume the non-Gaussianity is small and can be treated in perturbation theory. For large deviations from Gaussianity it is not trivial to construct a fully self-consistent theory. For example, if were locally power-law at each “step” in -space in a random walk, the resulting evaluated on each scale would no longer be a power-law; some violation of locality would be required so the distribution could “self-correct.” In any case, if we simply assume some pre-specified at all scales, it is still straightforward to evaluate the first-crossing distribution. The distribution in is given by the solution to the integro-differential equation:

 f(S) =−P′(δc|S)dδc(S)dS−∫δc(S)−∞∂P′(δ|S)∂Sdδ (30) P′(δ|S) ≡P(δ|S)−∫S0dS′f(S′)P(δ−δc[S′]|S−S′) (31)

where . This is essentially just the collapsed mass given by , corrected by the probability that the collapse occurred on a larger scale (smaller ), and can be solved numerically for any specified .

Consider the following form for the density PDF:

 dp(ρ)dlnρ∝exp(−γ|ln[ρ/¯ρ]|)={(ρ/¯ρ)γρ<¯ρ(ρ/¯ρ)−γρ≥¯ρ (32)

where . The exact functional form is arbitrary, of course, but convenient because it is a pure power-law symmetric in , and has a well-defined variance: We can therefore map this one-to-one to our assumed density power spectrum by assuming , with . Note that this gives over much of the dynamic range of interest, quite similar to the best-fit distributions in the references above. At low and high masses, the predicted mass function is slightly more shallow than our standard model. At high masses this is because of the more extended power-law tail to high ; at low masses this is both an effect of more first-crossings at larger scales and a result of some of the mass being moved from the “core” of the distribution to those tails. However, the differences are quite small. This is because a lognormal (unlike a pure normal distribution) is very similar to a single power law over a wide dynamic range. Moreover, the collapsed mass fraction is not extremely small, so it is not sampling some extreme tail of the distribution. So, for the same variance , deviations from lognormal behavior have only small effects.

### 3.5 The Core Mass Function

In this paper, we choose to focus on the mass function of GMCs and other large-scale structures in the ISM. Part of the reason for this is that we can focus on the first-crossing distribution (the largest scales on which structures are self-gravitating) and so have a well-defined mass function. Although there are certain similarities, this is not the same as the mass function of self-gravitating dense cores within GMCs, as calculated using qualitatively similar arguments in e.g. Padoan & Nordlund (2002) and Hennebelle & Chabrier (2008).

In principle, our model can be extended iteratively to smaller scales to investigate the mass function of cores and make a direct comparison with these previous predictions as well as observations, and in a companion paper (Hopkins, 2012) we attempt to do so. This is not trivial, however. The difficulty is that, because cores are substructures, the mass function definition (the resolution to the “cloud in cloud” problem) is somewhat ambiguous: we cannot simply isolate first-crossing. Even in simulations where the full three-dimensional properties are known, it is not trivial to find a unique mass function of such substructure in a turbulent medium (see e.g. Ballesteros-Paredes et al., 2006; Anathpindika, 2011). The approach of Hennebelle & Chabrier (2008) is to treat this ambiguity as an effective normalization term (and to truncate the problem at larger scales – treating the properties of the “parent” GMC as assumed/given and restricting to much smaller spatial scales); as such their derivation is similar to the original Press & Schechter (1974) derivation as discussed in § 1. That in Padoan & Nordlund 2002 more simply makes some general scaling arguments. But as we show in Fig. 1, this is not necessarily a good approximation. We therefore require some more detailed criteria to inform our definition of cores, for example some estimate of the scales on which fragmentation below the core scale will not occur (defining the “last-crossing,” as opposed to “first-crossing” distribution). This is a topic of considerable interest, but is outside the scope of our comparisons here.

## 4 Size-Mass & Linewidth-Size Relations

We can also use our model to predict the scaling laws obeyed by GMCs “at collapse.”

The linewidth-size relation follows trivially from our assumed turbulent power spectrum. The exact relation is plotted in Figure 3 for power-law turbulent slopes of and , with the normalization set by requiring a marginally stable disk with MW-like surface density. We can define the line width either as just the turbulent width, or the turbulent width plus the contribution from disk shear ; the distinction is unimportant for typical observed scales, but shear is predicted to contribute significantly to the velocities of the largest GMCs when . We compare with observations compiled from the MW and other local group galaxies from Bolatto et al. (2008) and Heyer et al. (2009).10

In the regime above the sonic length and below the scale height, this is just a simple power law with , i.e.  for or for . This is essentially an assumption of our model (although it follows from basic turbulent conditions); more interesting is that the normalization can be predicted from the assumption of marginal stability (), giving

 σv(R)≈0.4kms−1(⟨Σdisk⟩10M\sunpc−2)0.5(Rpc)0.5 (33)

This agrees well with the observed relation. In the full solution, because of the change in dimensionality above the scale , the relationship flattens if we consider only turbulent velocities; it becomes steeper, however, with the inclusion of the shear term.

This model also specifically predicts a residual dependence in the normalization of the linewidth-size relation that scales as , where is the large-scale mean disk surface density. We stress that this is not necessarily the same as a dependence on the local cloud (over a wide dynamic range, in fact, , hence , is similar). This is also, by definition, for bound objects, not for un-bound overdensities on small scales. The predicted dependence is shown indirectly in Figure 3, and directly in Figure 5, where we compare with the observations compiled in Heyer et al. (2009) in local galaxies and in Swinbank et al. (2011) for massive star-forming molecular complexes in lensed, high-redshift galaxies. These sample extremely different environments, and are indeed offset in the linewidth-size relation. However, the magnitude of their offsets is in good agreement with that predicted here.11 The galaxies in Swinbank et al. (2011) have an average surface density of , and a correspondingly very large measured (as expected for ); normalizing the predicted linewidth-size relation for these properties, we expect an order of magnitude larger at fixed size. Clouds observed in the MW center (Oka et al., 2001), which has a higher mean surface density than the local neighborhood but generally lower than estimated for the high-redshift systems, lie neatly between the predicted curves for the local and high-redshift cases (a mean offset of relative to the local clouds, corresponding to a factor of higher , about what is expected for the observed exponential profile). Similar offsets are known in other local galaxies with high surface densities, such as mergers and starburst galaxies (Wilson et al., 2003; Rosolowsky & Blitz, 2005).

As discussed in Hopkins et al. (2012), a dependence of exactly this sort is seen in high-resolution hydrodynamic simulations as well. In the observations, this normalization dependence has sometimes been interpreted as a consequence of magnetic support or confining external pressure (see the discussion in Blitz & Rosolowsky, 2006; Bolatto et al., 2008; Heyer et al., 2009), but in this context magnetic fields and pressure confinement are not explicitly present – such a scaling is a much more broad consequence of the simple Jeans requirements for collapse in any marginally stable environment.

The size-mass relation follows from the critical density derived in § 2.1, by simply inverting Eqn 14. We plot the exact prediction in Figure 4. In the regime above the sonic length but below the disk scale-height, recall that a power-law turbulent cascade gives the simple condition , so , i.e.  for , very similar to the observed power-law scaling. The normalization also follows – for MW-like global conditions

 Rcloud(R≫ℓsonic,R≪h)≈1.4σ−10.4pc(Mcloud300Msun)0.5 (34)

where is the normalization of the turbulent velocities . This corresponds to an approximately constant cloud surface density in agreement with Larson’s laws: in projection at the time of collapse. Note that re-calculating this for only changes the slope from to , well within the observational uncertainty. This will also alter behavior at the highest masses, but this is not significant until well above the mass function break. There does however appear to be tentative evidence for such a transition in the observations shown in Figure 4. As expected from the behavior of the linewidth-size relation, clouds in high density environments – which will have a higher in Eqn. 34 above – are offset to lower at fixed ; we show the same model prediction for the high-redshift systems in good agreement with the observations. Once again, MW center clouds and other local systems in environments with higher densities are similarly offset.

As discussed in § 3.5, fully extending the models here to the scales of dense cores is beyond the scope of this paper. However, we expect these cores, if self-gravitating, to obey the scaling in Figure 5. This means that if they form inside of high-density GMCs, we can (approximately) think of the “parent” GMC surface density as similar to the background term in Eqn. 33, and might expect them to have higher dispersions at fixed sizes. This has been seen in suggested from observations (Ballesteros-Paredes et al., 2011a), as part of a quasi-hierarchical gravitational collapse, similar to the predictions here. Of course, some regions can have much higher at fixed and be simply not self-gravitating; these will not lie on the relation in Figure 5 (they will be offset to higher ). This may, in turn, give rise to a dependence of the linewidth-size relation on the tracers and extinction threshold adopted, as observed (Goodman et al., 1998; Lombardi et al., 2010)

## 5 Spatial Clustering of GMCs

In analogy to dark matter halos, we can use the excursion set formalism to also determine the spatial clustering and correlation function strength of these bound sub-units. Following Mo & White (1996), the excess abundance of collapsing objects (relative to the mean abundance) in a sphere of radius with mean density is

 δcoll(R1,δc,1|R0,δ0)≡N(1|0)n(M1)V0−1 (35)

where is the average abundance of objects of mass (from the mass function) and is the number of collapsing objects in a region of radius (variance ) with fixed overdensity .

### 5.1 Linear Bias

If were constant, can be determined analytically and is simply

 N(1|0)=ρc,1V0M1δc,1−δ0√2π(S1−S0)3/2exp[−(δc,1−δ0)22(S1−S0)]dS1dM1 (36)

(Bond et al., 1991). In the regime where , so , this simplifies to

 (37)

where is defined as the linear bias of objects of mass .12

The barrier here is not constant. However, for arbitrary , we can calculate exactly by repeating our Monte Carlo excursion from § 3.1, but instead of beginning with initial conditions , for each walk, we begin at scale with density . The bias is then just the ratio of for small .

Figure 6 plots the bias as a function of cloud mass. A couple of key properties are clear. At high masses above the exponential cutoff in the mass function, the bias increases rapidly. This is qualitatively similar to what is seen for dark matter halos: because such systems are exponentially rare, they will tend to be strongly biased towards the few regions with substantial large-scale over densities. Physically, this corresponds to gas overdensities in the disk on scales larger than the scale-height , i.e. a preferential concentration of the most massive GMCs in global instabilities such as spiral arms, bars, and kpc-scale massive star-forming complexes, rather than their being randomly distributed across the gas. At intermediate masses below the mass function break, where most of the cloud mass lies, the bias is weak (order unity), so most of the mass in clouds simply traces most of the gas mass in general. We stress that this does not necessarily mean clouds are randomly distributed over the disk as a whole; it means they are unbiased relative to the gas mass distribution. But at low masses, the bias again rises (weakly). This is related to the anti-hierarchical nature of cloud formation: the bias here is driven by clouds which form via fragmentation from other clouds.

We can approximate these exact results using our previous approximate fitting functions for the mass function (Equation 2327) modified (as with the case of a linear barrier) so and . Neglecting the clouds-in-clouds problem (i.e. including those clouds), we obtain the approximate

 bcic(M1)≈ν21−1+ν1dlnρ1/dσ1δc,1(1+ν−11dlnρ1/dσ1) (38)

Which in practice is a small () correction to Equation 37. If we exclude clouds-in-clouds,

 b(M1)≈1δc,1[ν21−δc,1~B] (39)

(where is defined in Equation 23). This is identical to Eqn. 37 at high masses, but it allows for negative bias at low masses, if and . Physically, the fact that Equation 38 is always positive means that the number of bound regions of mass inside a large-scale overdensity always increases with . However, for some values of and , increasing more rapidly increases the probability that these regions are themselves inside a larger collapsed region. For a more detailed discussion of the leading-order corrections when considering a moving as opposed to constant barrier , we refer to Sheth et al. (2001).

### 5.2 The Correlation Function: Theory

Recall, the physical over-density is . The correlation function between collapsed objects of mass and background mass, as a function of radius , is defined by

 1+ξcm(R0,M1) ≡⟨N(1|R0)|ρ⟩n(M1)V0ρ0 (40) =⟨(1+δcoll)|exp(δ0−S[R0]/2)⟩R0 (41) =∫N(1|0)n(M1)V0e(δ0−S0/2)q(δ0|S0)dδ0 (42)

where the integral is over all , and is a weighting factor defined in Bond et al. (1991) as the probability that the overdensity at a random point, smoothed on a scale , is and does not exceed