THE STAR FORMATION RATE OF TURBULENT MAGNETIZED CLOUDS: COMPARING THEORY, SIMULATIONS, AND OBSERVATIONS

The PDF of the gas density in a turbulent medium—such as a molecular cloud—is the key ingredient for analytic models of star formation. A log-normal density PDF has been used to explain the mass distribution of cores and stars, the core mass function (CMF) and the stellar initial mass function (IMF; Padoan & Nordlund 2002; Hennebelle & Chabrier 2008, 2009; Elmegreen 2011; Veltchev et al. 2011; Donkov et al. 2012; Parravano et al. 2012; Hopkins 2012), the Kennicutt–Schmidt relation (Krumholz & McKee 2005; Tassis 2007), the SFE (Elmegreen 2008), and the SFR (Krumholz & McKee 2005; Padoan & Nordlund 2011; Hennebelle & Chabrier 2011). Here we concentrate on the SFR and derive its basic dependences.

The log-normal PDF of the gas density is defined as

$$\begin{equation} p_s(s)=\frac{1}{\sqrt{2\pi \sigma _s^2}}\exp \left(-\frac{(s-s_0)^2}{2\sigma _s^2}\right), \end{equation} \tag{ 1 }$$

expressed in terms of the logarithmic density,

$$\begin{equation} s\equiv \ln {(\rho /\rho _0)}. \end{equation} \tag{ 2 }$$

The PDF is a normal (Gaussian) distribution in s, meaning that it is a log-normal distribution in ρ. The quantities ρ₀ and s₀ denote the mean density and mean logarithmic density, the latter of which is related to the standard deviation σ_s by

$$\begin{equation} s_0=-\frac{1}{2}\,\sigma _s^2 \end{equation} \tag{ 3 }$$

due to the normalization and mass-conservation constraints of the PDF (Vázquez-Semadeni 1994; Federrath et al. 2008b). The reason to use s instead of ρ in the context of the density PDF is that s is dimensionless, and that the PDF of s is Gaussian unlike the PDF of ρ. This is because the distribution of ρ is generated by a multiplicative process in which shocks are amplified by other shocks as they collide and interact in isothermal supersonic turbulence, with the local Mach number being independent of the local density (Vázquez-Semadeni 1994; Passot & Vázquez-Semadeni 1998; Kritsuk et al. 2007; Federrath et al. 2010b). Since s∝ln (ρ) as defined in Equation (2), this multiplicative process in ρ turns into an additive process in s. Following the central limit theorem, a large sum of random variables produces a Gaussian distribution, and thus only p_s is Gaussian, while p_ρ is not. However, p_s can be easily transformed into p_ρ, because p_s ds = p_ρ dρ, and thus p_ρ = p_s/ρ (Li et al. 2003). We will omit the index s in p_s in the following and simply use p(s) for the PDF given by Equation (1).

As soon as significant collapse sets in, the density PDF develops a power-law tail at high densities (e.g., Klessen 2000; Kainulainen et al. 2009), which is discussed in more detail in Section 7.1.1 below and in Paper II (Federrath & Klessen 2012).

The standard deviation σ_s in Equation (1), which is a measure of how much the density varies in a turbulent medium, depends on (1) the amount of compression induced by the turbulent forcing mechanism, (2) the Mach number, and (3) the degree of magnetization. First, the turbulent energy injection mechanism determines the amount of compression induced directly by driving turbulence in the interstellar medium (ISM). Various turbulent driving mechanisms have been discussed and compared in Mac Low & Klessen (2004). For instance, expanding supernova shells (Balsara et al. 2004; de Avillez & Breitschwerdt 2005; Tamburro et al. 2009) or growing H ii regions around massive stars and clusters of stars (McKee 1989; Krumholz et al. 2006; Gritschneder et al. 2009; Peters et al. 2010; Goldbaum et al. 2011) as well as compression of ISM gas in galactic spiral shocks (Elmegreen 2009) and gravitational contraction (Hoyle 1953; Vazquez-Semadeni et al. 1998; Klessen & Hennebelle 2010; Elmegreen & Burkert 2010; Federrath et al. 2011c) are likely exciting a considerable amount of compressible modes that will directly lead to compression, and thus to higher density contrasts on molecular cloud scales in the ISM, while galactic rotation and magnetorotational instabilities (e.g., Piontek & Ostriker 2004, 2007) are likely producing more solenoidal modes. Second, higher Mach numbers $\mathcal {M}$ lead to stronger shocks and thus to higher density contrasts. For instance, the density jump in a non-magnetized, isothermal shock is proportional to $\mathcal {M}^2$. Finally, higher magnetization dampens density fluctuations as magnetic fields act like a cushion due to the additional magnetic pressure (Ostriker et al. 2001; Price et al. 2011).

The actual dependence of turbulent density fluctuations σ_s on the three parameters above (forcing, Mach number, and magnetic field) can be derived from the shock jump conditions of an individual MHD shock, and then averaged over a whole ensemble of such shocks (Padoan & Nordlund 2011). Molina et al. (2012) provide a rigorous derivation of σ_s for different correlations of the magnetic field with density. They distinguish three cases, B∝ρ⁰, B∝ρ^1/2, and B∝ρ¹. For the intermediate case, Molina et al. (2012) derive

$$\begin{equation} \sigma _s^2=\ln \left(1+b^2\mathcal {M}^2\frac{\beta }{\beta +1}\right)\,, \end{equation} \tag{ 4 }$$

which is similar to the relation derived in Padoan & Nordlund (2011), except for the factor b², explained below, and except for the definition of β, for which Padoan & Nordlund (2011) only take post-shock gas into account (see the more extended discussion on this issue in Section 2.4.2). The case B∝ρ¹ is similar to the intermediate case, but is a rather extreme MHD case because magnetic-field lines are assumed to be oriented only perpendicular to the flow direction; so is the other extreme case in which the magnetic field is assumed to be parallel to the flow, yielding B∝ρ⁰. In the more realistic case of turbulent flows, field lines become tangled, and the B–ρ correlation is a combination of compression of field lines and turbulent dynamo amplification (Schleicher et al. 2010; Sur et al. 2010; Federrath et al. 2011c; Turk et al. 2012; Schober et al. 2012). In a three-dimensional system with a random distribution of flow velocities and magnetic-field orientations, B∝ρ^1/2 provides a reasonable intermediate dependence. We will thus only consider B∝ρ^1/2 here, which is favored by simulations (Padoan & Nordlund 1999; Collins et al. 2011; Molina et al. 2012) and is also close to what is suggested from observations of magnetic fields in molecular clouds (Crutcher et al. 2010).³

In the case of B∝ρ⁰, i.e., for no density correlation of the magnetic field, Equation (4) reduces to the well-known and frequently used hydrodynamic (HD) expression, $\sigma _s^2=\ln \left(1+b^2\mathcal {M}^2\right)$ with β → ∞ (e.g., Padoan et al. 1997; Passot & Vázquez-Semadeni 1998; Ostriker et al. 2001; Lemaster & Stone 2008; Federrath et al. 2008b; Price et al. 2011) as a necessary condition in the purely HD limit. The parameters b, $\mathcal {M}$, and β in Equation (4) are the turbulent forcing parameter, the rms sonic Mach number, and the ratio of thermal-to-magnetic pressure, plasma β = P_th/P_mag. Using the definitions of the thermal pressure for an isothermal equation of state P_th = ρc²_s, magnetic pressure P_mag = B²/(8π), Alfvén velocity v²_A = B²/(4πρ), sonic and Alfvén Mach numbers, $\mathcal {M}=\sigma _V/c_{\rm s}$ and $\mathcal {M}_{\rm A}=\sigma _V/v_{\rm A}$, the plasma beta can be expressed as $\beta =2c_{\rm s}^2/v_{\rm A}^2=2\mathcal {M}_{\rm A}^2/\mathcal {M}^2$. These are all dimensionless numbers, rendering them particularly useful because they determine the basic properties of turbulent plasmas and can thus be directly compared for any such system. Equation (4) can thus also be written as

$$\begin{equation} \sigma _s^2=\ln \left(1+b^2\mathcal {M}^2\frac{2\mathcal {M}_{\rm A}^2}{\mathcal {M}^2+2\mathcal {M}_{\rm A}^2}\right)\,. \end{equation} \tag{ 5 }$$

The forcing parameter b was shown to vary smoothly between b ≈ 1/3 for purely solenoidal (divergence-free) forcing and b ≈ 1 for purely compressive (curly-free) forcing of the turbulence (Federrath et al. 2008b, 2010b; Schmidt et al. 2009; Seifried et al. 2011b; Micic et al. 2012; Konstandin et al. 2012a). A stochastic mixture of forcing modes in three-dimensional space leads to b ≈ 0.4 (see Figure 8 in Federrath et al. 2010b).

Using numerical simulations, Molina et al. (2012) found that Equations (4) and (5) work well in the regime $\mathcal {M}_{\rm A}\gtrsim 2$, while for $\mathcal {M}_{\rm A}\lesssim 2$, the assumption of isotropy entering the analytic derivation of Equations (4) and (5) breaks down, so we only apply them in the super-Alfvénic regime throughout the following.

Here we present an analytic derivation of the SFR from the statistics of supersonic, isothermal magnetized turbulence. The main ingredient for this analytic derivation is an integral over the density PDF, Equation (1), in order to estimate the gas mass above a given density threshold, contributing to star formation. We will compare different ways of estimating the density threshold, which is the main difference between the three most successful existing analytic models for the SFR (Krumholz & McKee 2005; Padoan & Nordlund 2011; Hennebelle & Chabrier 2011). We will express all quantities in terms of dimensionless numbers, in order to simplify the derivation and to make it more general. We follow the standard terminology and use the star formation rate per freefall time (SFR_ff), as coined by Krumholz & McKee (2005), which is the mass fraction accreting onto stars per unit time, where the time is expressed in units of the mean freefall time.

The SFR in units of M_☉ yr⁻¹ can be computed by scaling SFR_ff with the real cloud mass M_c and the actual freefall time evaluated at the mean density of the cloud, t_ff(ρ₀):

$$\begin{equation} \mathrm{SFR} \equiv \frac{M_{\rm c}}{t_{\rm ff}(\rho _0)}\,\mathrm{SFR}_{\rm ff}\,. \end{equation} \tag{ 6 }$$

Note that this definition of SFR_ff is different from the definition used in Krumholz & Tan (2007) and Krumholz et al. (2012), who use freefall times estimated at different densities and/or use a definition based on column densities, such that the values of SFR_ff quoted in those studies and the ones computed here cannot be directly compared. For instance, given an SFR for fixed M_c, the dimensionless value of SFR_ff would be much smaller, if the freefall time at a high-density tracer was used rather than the freefall time at the mean density of the cloud, because t_ff(ρ > ρ₀) is shorter than t_ff(ρ₀).

The basic idea for an analytic model of SFR_ff is to integrate the log-normal density PDF, Equation (1), weighted by ρ/ρ₀ to get the mass fraction of gas with density above a critical density s_crit (to be determined below in Section 2.4) and weighted by a freefall-time factor to construct a dimensionless mass rate:

$$\begin{equation} \mathrm{SFR}_{\rm ff}= \frac{\epsilon }{\phi _t} \int _{s_{\rm crit}}^{\infty }{\frac{t_{\rm ff}(\rho _0)}{t_{\rm ff}(\rho)} \frac{\rho }{\rho _0} \, p(s) \,{d}s}\,. \end{equation} \tag{ 7 }$$

Note that the factor t_ff(ρ₀)/t_ff(ρ) appears inside the integral because gas with different densities has different freefall times,

$$\begin{equation} t_{\rm ff}(\rho)\equiv \left(\frac{3\pi }{32G\rho }\right)^{1/2}\,, \end{equation} \tag{ 8 }$$

which should be taken into account in the most general case (see Hennebelle & Chabrier 2011). Previous estimates for SFR_ff either used a factor t_ff(ρ₀)/t_ff(ρ₀) = 1 (Krumholz & McKee 2005) or a factor t_ff(ρ₀)/t_ff(ρ_crit) with ρ_crit = ρ₀exp (s_crit) (Padoan & Nordlund 2011), both of which are independent of density and were thus pulled out of the integral. We will show, however, that it is crucial to take the multi-freefall nature of gas with different densities into account to obtain better models for SFR_ff.

The constant factor in Equation (7) accounts for the fact that only a certain fraction of the gas above s_crit might actually accrete onto stars. Since individual stars form in accretion disks from which powerful jets, winds, and outflows are launched during the process of stellar birth, it is likely that a certain fraction of the accreted material is re-injected into the ISM, thus leading to < 1. Theoretical upper limits are in the range ≈ 0.25–0.7 (e.g., Matzner & McKee 2000). The observed displacement of the characteristic mass in the IMF (e.g., Kroupa 2001; Chabrier 2003) with respect to the CMF (e.g., Johnstone et al. 2000) has been taken to argue that might be around 0.3–0.5 (Alves et al. 2007; André et al. 2010); see however Ward et al. (2012).

The factor 1/ϕ_t in Equation (7) is also of the order of unity and accounts for the uncertainty in the timescale factor t_ff(ρ₀)/t_ff(ρ), originally introduced in Krumholz & McKee (2005). We will determine the best-fit values of and 1/ϕ_t in Sections 4 and 6, when we compare the theories with simulations and observations.

In the following, we will solve Equation (7), using different density thresholds s_crit, according to the previous analytic studies of the SFR by Krumholz & McKee (2005, KM), Padoan & Nordlund (2011, PN), and Hennebelle & Chabrier (2011, HC).⁴ We distinguish six cases, named "KM," "PN," "HC," and "multi-ff KM," "multi-ff PN," "multi-ff HC," as distinguished in Hennebelle & Chabrier (2011). The first three represent the original analytic derivations by Krumholz & McKee (2005), Padoan & Nordlund (2011), and Hennebelle & Chabrier (2011), while the last set of three are all based on the multi-freefall expression of the integral (7). The difference for this last set of three is only the model for the critical density, i.e., the lower limit of the integral. We note that the ideas inherent in each of the original theories contribute to our present understanding of the turbulence-regulated SFR. Krumholz & McKee (2005) developed the basic model, Padoan & Nordlund (2011) extended it to include magnetic fields, and Hennebelle & Chabrier (2011) improved all models by introducing multi-freefall versions of the aforementioned theories, yet without considering magnetic fields. We build on all these approaches and extend the non-magnetic multi-freefall models to include magnetic fields. We then determine the best combination of the aforementioned theoretical ideas to come up with a more universal theoretical model for the SFR. Table 1 summarizes all six theoretical models, which are discussed and derived in detail in the following.

Table 1. Six Analytic Models for the Star Formation Rate per Freefall Time

Analytic Model	Freefall-time Factor	Critical Density ρcrit/ρ0 = exp (scrit)	SFRff
KM	1	(π2/5) ϕ2x × $\alpha _{\rm vir}\mathcal {M}^2(1+\beta ^{-1})^{-1}$	/(2ϕt){1 + erf[(σ2s − 2scrit)/(8σ2s)1/2]}
PN	tff(ρ0)/tff(ρcrit)	(0.067) θ−2 × $\alpha _{\rm vir}\mathcal {M}^2 f(\beta)$	/(2ϕt){1 + erf[(σ2s − 2scrit)/(8σ2s)1/2]}exp [(1/2)scrit]
HC	tff(ρ0)/tff(ρ)	(π2/5) y−2cut × $\alpha _{\rm vir}\mathcal {M}^{-2}(1+\beta ^{-1}) + \tilde{\rho }_{\rm crit,turb}$	/(2ϕt){1 + erf[(σ2s − scrit)/(2σ2s)1/2]}exp [(3/8)σ2s]
multi-ff KM	tff(ρ0)/tff(ρ)	(π2/5) ϕ2x × $\alpha _{\rm vir}\mathcal {M}^2(1+\beta ^{-1})^{-1}$	/(2ϕt){1 + erf[(σ2s − scrit)/(2σ2s)1/2]}exp [(3/8)σ2s]
multi-ff PN	tff(ρ0)/tff(ρ)	(0.067) θ−2 × $\alpha _{\rm vir}\mathcal {M}^2 f(\beta)$	/(2ϕt){1 + erf[(σ2s − scrit)/(2σ2s)1/2]}exp [(3/8)σ2s]
multi-ff HC	tff(ρ0)/tff(ρ)	(π2/5) y−2cut × $\alpha _{\rm vir}\mathcal {M}^{-2}(1+\beta ^{-1})$	/(2ϕt){1 + erf[(σ2s − scrit)/(2σ2s)1/2]}exp [(3/8)σ2s]

Notes. The function f(β), entering the critical density in the PN and multi-ff PN models, is given by Equation (31). The added turbulent contribution $\tilde{\rho }_{\rm crit,turb}$ in the critical density of the HC model is given by Equation (39).

Download table as:  ASCIITypeset image

2.4.1. The KM Model

In the KM model by Krumholz & McKee (2005), the freefall-time factor t_ff(ρ₀)/t_ff(ρ) in Equation (7) is simply set to unity. Moreover, Krumholz & McKee (2005) define the critical density s_crit in the lower limit of the SFR_ff integral by comparing the Jeans (1902) length

$$\begin{equation} \lambda _{\rm J}(\rho)=\left(\frac{\pi c_{\rm s}^2}{G\rho }\right)^{1/2}, \end{equation} \tag{ 9 }$$

evaluated at the mean density with the sonic scale λ_s (defined in Equation (13) below),

$$\begin{equation} s_{\rm crit}=2\,\ln \left(\phi _x\frac{\lambda _{\rm J}(\rho _0)}{\lambda _{\rm s}}\right). \end{equation} \tag{ 10 }$$

This choice is motivated by the expectation that the collapse sets in roughly at the sonic scale, where the turbulent fluctuations are of the order of the thermal sound speed, i.e., the local Mach number has dropped to about unity at the sonic scale (Vázquez-Semadeni et al. 2003; Federrath et al. 2010b). The global turbulent supersonic support is expected to become insignificant at the sonic scale, such that collapse can proceed below that scale (e.g., Mac Low & Klessen 2004). The leading factor 2 in Equation (10) stems from the density dependence of the Jeans length, λ_J(ρ)∝ρ^−1/2, and the numerical factor ϕ_x allows for slight variations in the actual scale on which the collapse sets in. Krumholz & McKee (2005) find ϕ_x = 1.12 for the simulations by Vázquez-Semadeni et al. (2003). In real molecular clouds, the sonic scale is expected to be of the order of 0.1 pc within factors of a few (e.g., Falgarone et al. 1992; Goodman et al. 1998; Stahler & Palla 2004; Schnee et al. 2007; McKee & Ostriker 2007).

To make Equation (10) more useful, we express all dependent variables for s_crit in terms of dimensionless numbers. This can be achieved by rewriting the Jeans length as

$$\begin{equation} \lambda _{\rm J}(\rho _0)=\left(\frac{\pi c_{\rm s}^2}{G\rho _0}\right)^{1/2}=\pi c_{\rm s}\left(\frac{L^3}{6GM_{\rm c}}\right)^{1/2}\,, \end{equation} \tag{ 11 }$$

where we have assumed a spherical cloud with diameter L, mass M_c, and isothermal sound speed c_s. Since the velocity fluctuations in a turbulent medium depend on the length scale ℓ as

$$\begin{equation} \sigma _v(\ell)=\sigma _V\,(\ell /L)^p \,, \end{equation} \tag{ 12 }$$

where σ_V ≈ 1 km s⁻¹ is the three-dimensional, non-thermal velocity dispersion on the scale L ≈ 1 pc, and p ≈ 0.5 from observations in Galactic clouds (Larson 1981; Solomon et al. 1987; Ossenkopf & Mac Low 2002; Heyer & Brunt 2004; Heyer et al. 2009; Roman-Duval et al. 2011), the Galactic Central Molecular Zone (Jones et al. 2012; Shetty et al. 2012), and from numerical simulations (Kritsuk et al. 2007; Schmidt et al. 2009; Federrath et al. 2010b), the sonic scale can be written as

$$\begin{equation} \lambda _{\rm s}=L\left(c_{\rm s}/\sigma _V\right)^{1/p}\,. \end{equation} \tag{ 13 }$$

Substituting Equations (11) and (13) into Equation (10), we find

$$\begin{eqnarray} s_{\rm crit}{_{\rm ,KM}} \,& = &\, \ln \left[\frac{\phi _x^2\pi ^2}{5}\,\frac{5\sigma _V^2 L}{6{\it GM}_{\rm c}}\,\left(\frac{\sigma _V}{c_{\rm s}}\right)^{2(1-p)/p}\right] \nonumber \\ \,& = &\, \ln \left[(\pi ^2/5)\phi _x^2\,\alpha _{\rm vir}\,\mathcal {M}^{2}\right]\,, \end{eqnarray} \tag{ 14 }$$

where we have identified the virial parameter for a spherical, uniform-density cloud with velocity dispersion σ_V on the diameter scale L,

$$\begin{equation} \alpha _{\rm vir,\circ }=5\sigma _V^2 L/(6GM_{\rm c}), \end{equation} \tag{ 15 }$$

and the rms Mach number, $\mathcal {M}=\sigma _V/c_{\rm s}$, and used p = 0.5 in the second step. This derivation is essentially identical to the one presented in Krumholz & McKee (2005), with the exception that we use the more general expression for the virial parameter here,

$$\begin{equation} \alpha _{\rm vir}= 2E_{\rm kin}/|E_{\rm grav}|, \end{equation} \tag{ 16 }$$

the ratio of twice the kinetic energy to the gravitational energy. This general form reduces to α_vir, ° given by Equation (15) with E_kin = M_cσ²_V/2 and E_grav = −3GM²_c/(5R) for a spherical, homogeneous cloud with radius R = L/2. We emphasize that the definition of α_vir, ° is based on global parameters, assuming a spherical cloud with uniform density. This is far from realistic, given that clouds are in fact highly inhomogeneous and non-spherical. In the analytic derivations, however, this simplification given by Equation (15) is necessary to enable a mathematical analysis of the problem. In the simulations discussed in Section 3 below, however, we will directly compute the virial parameter from the gravitational potential of the actual, three-dimensional, inhomogeneous spatial gas distribution, providing a more general and accurate measure of the virial parameter given by the general form, Equation (16). This is discussed further below when we compare the theories to numerical simulations and in Section 7.1.3.

The original model by Krumholz & McKee (2005) neglects magnetic fields. Here, magnetic-field effects are partially added automatically by using Equation (4) for σ_s, such that σ_s decreases with increasing magnetic energy, as derived in Molina et al. (2012). This, however, only changes σ_s, while a modification of s_crit is also necessary to fully account for magnetic-pressure effects on SFR_ff.

Here we provide and apply a simple rule to include magnetic-field effects in the expression for the critical density s_crit. The key idea is to replace the thermal pressure by the sum of the thermal and magnetic pressures:

$$\begin{eqnarray} P_{\rm th} & \rightarrow & P_{\rm th} + P_{\rm mag} \nonumber \\ \iff \rho c_{\rm s}^2 & \rightarrow & \rho c_{\rm s}^2 + (1/2)\rho v_{\rm A}^2, \end{eqnarray} \tag{ 17 }$$

where the second line implies isothermal gas. Using v²_A = 2c_s²β⁻¹ with the definition of plasma β = P_th/P_mag in Section 2.2, we can thus simply replace the sound speed by an effective sound speed,

$$\begin{equation} c_{\rm s}\rightarrow c_{\rm s}(1+\beta ^{-1})^{1/2}. \end{equation} \tag{ 18 }$$

Since $\mathcal {M}=\sigma _V/c_{\rm s}$, we can also replace the sonic Mach number by an effective Mach number to take magnetic pressure into account:

$$\begin{equation} \mathcal {M}\rightarrow \mathcal {M}(1+\beta ^{-1})^{-1/2}. \end{equation} \tag{ 19 }$$

Doing this for s_{crit, KM} in Equation (14) yields the magnetic version of the critical density,

$$\begin{equation} s_{\rm crit}{_{\rm ,KM}} = \ln \big[(\pi ^2/5)\phi _x^2\,\alpha _{\rm vir}\,\mathcal {M}^{2}(1+\beta ^{-1})^{-1}\big]. \end{equation} \tag{ 20 }$$

Even though we simply replaced the thermal sound speed by an effective, magnetic sound speed to derive this expression, it has a deeper physical meaning. What we physically do in the derivation of s_crit is to replace the thermal Jeans length in the numerator of Equation (10) with the magnetothermal Jeans length,

$$\begin{equation} \lambda _{\rm J,mag}= \left(\frac{\pi c_{\rm s}^2(1+\beta ^{-1})}{G\rho }\right)^{1/2}, \end{equation} \tag{ 21 }$$

and the sonic scale in the denominator with the magnetosonic scale,

$$\begin{equation} \lambda _{\rm ms}= L[c_{\rm s}(1+\beta ^{-1})^{1/2}/\sigma _V]^{1/p}. \end{equation} \tag{ 22 }$$

We note that the magnetic modifications given by Equations (17) only account for magnetic pressure, i.e., isotropic pressure induced by the small-scale magnetic field. It does not account for mean magnetic-field effects, and as such will only be a valid extension to MHD as long as the turbulence remains trans- to super-Alfvénic because sub-Alfvénic turbulence with a strong mean magnetic-field component is anisotropic, which is discussed at more detail below.

Finally, solving the general SFR_ff-integral (Equation (7)) with s_crit = s_{crit, KM} from Equation (20) and unity for the freefall-time factor (see Table 1 for a summary), the SFR per freefall time in the KM model is

$$\begin{eqnarray} \mathrm{SFR}_{\rm ff}{_{\rm ,KM}} \,& = &\, \frac{\epsilon }{\phi _t} \int_{s_{\rm crit}{_{\rm ,KM}}}^{\infty } {\exp (s)\,p(s) \,{d}s} \nonumber \\ \,& = &\, \frac{\epsilon }{2\phi _t} \left[1+\mathrm{erf}\left(\frac{\sigma _s^2-2s_{\rm crit}{_{\rm ,KM}}}{\sqrt{8\sigma _s^2}}\right)\right]. \end{eqnarray} \tag{ 23 }$$

This derivation is identical to the one in Krumholz & McKee (2005), except for the extension to include magnetic fields in the theory based on the plasma β terms in σ_s, Equation (4), and in the critical density, Equation (20).

2.4.2. The PN Model

Padoan & Nordlund (2011) use t_ff(ρ₀)/t_ff(ρ_crit) as the freefall-time factor t_ff(ρ₀)/t_ff(ρ) in Equation (7), such that the freefall time of the critical density is used for all densities above the critical density to estimate SFR_ff. Unlike Krumholz & McKee (2005) who relate the critical density s_crit to the Jeans length and the sonic scale, Padoan & Nordlund (2011) related the critical density to the magnetic shock jump conditions and to the magnetic critical mass for collapse. Starting with their assumed balance of thermal plus magnetic pressure by turbulent ram pressure on the cloud scale,

$$\begin{equation} \rho _{\rm MHD}\left(c_{\rm s}^2+\frac{1}{2}v_{\rm A}^2\right) = \rho _0\left(\frac{\sigma _V}{2}\right)^2, \end{equation} \tag{ 24 }$$

and using the definitions for $\mathcal {M}$ and β from Section 2.2, Padoan & Nordlund (2011) arrive at an expression for the density jump

$$\begin{equation} \rho _{\rm MHD} = \rho _0\frac{\mathcal {M}^2}{4}\frac{\beta }{\beta +1}. \end{equation} \tag{ 25 }$$

This leads to the post-shock thickness

$$\begin{equation} \lambda _{\rm MHD} = \theta L \, \frac{4}{\mathcal {M}^2}\frac{\beta +1}{\beta }, \end{equation} \tag{ 26 }$$

since ρ_MHD/ρ₀ = θL/λ_MHD with the numerical parameter θ ≲ 1, the fraction of the cloud size forming the largest shocks. Thus, θL can be interpreted as the turbulent injection or forcing scale. In numerical simulations, most of the kinetic, turbulent energy is usually injected at a wavenumber k = 2 in units of 2π/L, corresponding to half of the total cloud size (e.g., Kritsuk et al. 2007; Schmidt et al. 2009; Federrath et al. 2010b), as in the simulations discussed below in Section 3. Thus, θ ≈ 1/2, but there might be some corrections to that particular scale (Wang & George 2002). Padoan & Nordlund (2011) take θ ≈ 0.35. Here, we will simply interpret θ as a numerical factor of the order of unity, accounting for any uncertainty in the post-shock thickness with respect to the total cloud scale L in Equation (26).

In order to derive a critical density for star formation, Padoan & Nordlund (2011) compare the mass of a sphere with radius λ_MHD/2 to the critical mass for collapse. McKee (1989) define the critical mass for collapse of a magnetized gas sphere as

$$\begin{equation} M_{\rm crit}\approx M_{\rm BE}+ M_\Phi , \end{equation} \tag{ 27 }$$

where

$$\begin{equation} M_{\rm BE}= 1.182 c_{\rm s}^3 G^{-3/2} \rho ^{-1/2} \end{equation} \tag{ 28 }$$

is the Bonnor–Ebert mass (Ebert 1955; Bonnor 1956) and

$$\begin{equation} M_\Phi = m_\Phi \frac{\pi R^2\,B}{G^{1/2}} = m_\Phi ^3\frac{9\pi ^{5/2}}{2G^{3/2}}\rho ^{-1/2}v_{\rm A}^3 \end{equation} \tag{ 29 }$$

is the magnetic critical mass for a sphere with radius R, threaded by a magnetic field B, where we have used the Alfvén velocity v_A = B/(4πρ)^1/2 in the second step. The numerical factor $m_\Phi$ in Equation (29) can vary depending on the geometry and model taken, e.g., Padoan & Nordlund (2011) take $m_\Phi =0.17$ with a reference to Tomisaka et al. (1988), while McKee (1989) use $m_\Phi =0.12$, and Strittmatter (1966) derive $m_\Phi =(12\pi ^2/5)^{-1/2}\approx 0.21$ for a non-rotating cloud and $m_\Phi =(9\pi ^4/10)^{-1/2}\approx 0.11$ for an oblate spheroidal cloud with eccentricity approaching unity (see Nakano & Nakamura 1978).

Finally, inserting Equations (28) and (29) into Equation (27) and setting the critical mass M_crit(ρ_crit) = (4π/3)(λ_MHD/2)³ρ_crit with the post-shock thickness given by Equation (26) yields the critical density,

$$\begin{equation} s_{\rm crit}{_{\rm ,PN}} = \ln {[0.067\theta ^{-2} \alpha _{\rm vir}\mathcal {M}^2 f(\beta) ]} \end{equation} \tag{ 30 }$$

with

$$\begin{equation} f(\beta) \equiv \frac{(1+0.925 \beta ^{-3/2})^{2/3}}{(1+\beta ^{-1})^2}. \end{equation} \tag{ 31 }$$

Note that s_{crit, PN} has the same dependence on α_vir and $\mathcal {M}$ as s_{crit, KM} in Equation (20).

Padoan & Nordlund (2011) use a rather special definition of β, which is the average post-shock β. From a semi-analytical comparison of the mean magnetic field with the rms magnetic field, they derive a criterion for β based on the average Alfvén Mach number, which Padoan & Nordlund (2011) simply use as a switch between MHD and purely HD turbulence. However, it is not straightforward to derive a post-shock value of β because it involves a density-threshold dependence (see discussion in Padoan & Nordlund 2011). Moreover, the switch discussed by Padoan & Nordlund (2011) is a semi-analytical criterion derived from their simulations. We therefore decide to ignore this special definition of β for simplicity and apply Equation (30) with our definition of β (see Section 2.2), which includes all, and not just the post-shock, gas. This is consistent with the definition of all other dynamical quantities of interest, e.g., α_vir, $\mathcal {M}$, $\mathcal {M}_{\rm A}$, ρ₀, etc.

Using t_ff(ρ₀)/t_ff(ρ_{crit, PN}) and inserting s_{crit, PN} into the general Equation (7) for SFR_ff yields

$$\begin{eqnarray} &&\mathrm{SFR}_{\rm ff}{_{\rm ,PN}} = \frac{\epsilon }{\phi _t} \exp \left(\frac{1}{2}s_{\rm crit}{_{\rm ,PN}}\right) \int _{s_{\rm crit}{_{\rm ,PN}}}^{\infty } {\exp (s)\,p(s) \,{d}s} \nonumber \\ &&= \frac{\epsilon }{2\phi _t} \exp \left(\frac{1}{2}s_{\rm crit}{_{\rm ,PN}}\right) \left[1+\mathrm{erf}\left(\frac{\sigma _s^2-2s_{\rm crit}{_{\rm ,PN}}}{\sqrt{8\sigma _s^2}}\right)\right] \nonumber \\ \end{eqnarray} \tag{ 32 }$$

for the PN model.

2.4.3. The HC Model

Hennebelle & Chabrier (2011) were the first to argue that the freefall-time factor t_ff(ρ₀)/t_ff(ρ) must be used in Equation (7), such that different densities contribute to SFR_ff with their individual freefall time (see Equation (8)). The full HC model for SFR_ff is based on the mass spectrum of gravitationally bound structures, as derived in Hennebelle & Chabrier (2008, 2009):

$$\begin{equation} \mathcal {N}(M) = \frac{\,{d}(N/V)}{\,{d}M} \propto \,-\frac{1}{M}\frac{\,{d}s}{\,{d}M}\,\exp (s)\,p(s)\,, \end{equation} \tag{ 33 }$$

which is essentially Equation (6) in Hennebelle & Chabrier (2011), except for the freefall-time factor. The SFR in the HC model is then given by the integral over the mass spectrum, weighted by the mass and the freefall-time factor:

$$\begin{eqnarray} \mathrm{SFR}_{\rm ff}\,& = &\, -\frac{\epsilon }{\phi _t} \int_{0}^{M_{\rm cut}} \frac{M\,{d}M}{M}\frac{\,{d}s}{\,{d}M} \frac{t_{\rm ff}(\rho _0)}{t_{\rm ff}(\rho)}\,\exp (s)\,p(s) \nonumber \\ \,& = &\, \frac{\epsilon }{\phi _t} \int_{s_{\rm crit}}^{\infty }{\frac{t_{\rm ff}(\rho _0)}{t_{\rm ff}(\rho)} \frac{\rho }{\rho _0} \, p(s) \,{d}s}\,. \end{eqnarray} \tag{ 34 }$$

Note that the first equality is the same as Equation (7) in Hennebelle & Chabrier (2011).⁵ It can be simplified to the second line in Equation (34), by transforming the mass variable into the logarithmic density variable s and changing the limits of the integral accordingly. We emphasize that the second equality in Equation (34) is exactly the same as the general model for SFR_ff given by Equation (7) above.

In the HC model, the critical density s_{crit, HC} is defined by requiring that the turbulent Jeans length λ_{J, turb} at the critical density is a fraction y_cut of the cloud scale L. Hennebelle & Chabrier (2011) do not provide an explicit physical interpretation of this choice, but a follow-up study is in preparation (P. Hennebelle & G. Chabrier 2012, private communication). The turbulent Jeans length is obtained by adding an effective turbulent pressure (see Chandrasekhar 1951a, 1951b; Bonazzola et al. 1987)⁶ to the sound speed in the purely thermal Jeans length, Equation (9):

$$\begin{eqnarray} \lambda _{\rm J,turb}\,& \equiv &\, \left(\frac{\pi c_{\rm s}^2+(\pi /3)\sigma _v^2(\lambda _{\rm J,turb})}{G\rho }\right)^{1/2} \nonumber \\ \,& = &\, \left(\frac{\pi c_{\rm s}^2+\pi \lambda _{\rm J,turb}\sigma _V^2/(3L)}{G\rho }\right)^{1/2} \,, \end{eqnarray} \tag{ 35 }$$

in which the turbulent velocity dispersion, Equation (12), must be evaluated on the scale ℓ = λ_{J, turb}, such that the turbulent Jeans length is implicitly defined by Equation (35). Rewriting yields a quadratic equation with two solutions:

$$\begin{equation} \lambda _{\rm J,turb}(\rho)=\frac{\pi \sigma _V^2 \pm \sqrt{36\pi c_{\rm s}^2GL^2\rho +\pi ^2\sigma _V^4}}{6GL\rho }\, \end{equation} \tag{ 36 }$$

for which only the positive root is physical because the Jeans length must become larger when adding a stabilizing pressure—in this case a turbulent pressure. Naturally, this expression reduces to the thermal Jeans length for σ_V → 0. Now, setting the turbulent Jeans length equal to y_cutL as defined in Hennebelle & Chabrier (2011), and identifying the virial parameter, Equation (15), and the Mach number $\mathcal {M}=\sigma _V/c_{\rm s}$ finally yields the critical density threshold in the HC model:

$$\begin{equation} s_{\rm crit}{_{\rm ,HC}} = \ln {[ \tilde{\rho }_{\rm crit,th} + \tilde{\rho }_{\rm crit,turb} ]} , \end{equation} \tag{ 37 }$$

where the (magneto)thermal contribution is

$$\begin{equation} \tilde{\rho }_{\rm crit,th} \equiv (\pi ^2/5)y_{\rm cut}^{-2} \alpha _{\rm vir}\mathcal {M}^{-2} (1+\beta ^{-1}) \,, \end{equation} \tag{ 38 }$$

and the turbulent contribution is

$$\begin{equation} \tilde{\rho }_{\rm crit,turb} \equiv (\pi ^2/15)\,y_{\rm cut}^{-1}\,\alpha _{\rm vir}\,. \end{equation} \tag{ 39 }$$

Note that the dependence of the thermal contribution to s_{crit, HC} on α_vir is the same as in the KM and PN models, while the dependence on the Mach number is $\mathcal {M}^{-2}$, which is the opposite of the dependence in the KM and PN models, for both of which $\rho _{\rm crit}\propto \mathcal {M}^{+2}$ (see Table 1 for a summary of all analytic models).

The original HC model does not take magnetic fields into account, but we have extended the HC theory to MHD here by replacing the sonic Mach number in Equation (38) with the magnetic version in the same way as done for the KM model via Equation (19). The magnetic correction factor (1 + β⁻¹) in Equation (38) simply becomes unity in the hydrodynamical limit (β → ∞).

The SFR in the HC model is thus given by integrating Equation (34) or equivalently Equation (7) with s_{crit, HC}, which yields

$$\begin{eqnarray} \mathrm{SFR}_{\rm ff}{_{\rm ,HC}} & = & \frac{\epsilon }{\phi _t} \int_{s_{\rm crit}{_{\rm ,HC}}}^{\infty } {\exp \left(\frac{3}{2}s\right) \, p(s) \,{d}s} \nonumber \\ & = & \frac{\epsilon }{2\phi _t} \exp \left(\frac{3}{8}\sigma _s^2\right) \left[1+\mathrm{erf}\left(\frac{\sigma _s^2-s_{\rm crit}{_{\rm ,HC}}}{\sqrt{2\sigma _s^2}}\right)\right]. \nonumber \\ \end{eqnarray} \tag{ 40 }$$

2.4.4. The Multi-freefall KM Model

Following Hennebelle & Chabrier (2011), we define all three multi-freefall versions of the KM, PN, and HC models by solving the generalized, multi-freefall integral, Equation (7). The analytic solution of that equation for an arbitrary threshold s_crit is

$$\begin{equation} \mathrm{SFR}_{\rm ff}= \frac{\epsilon }{2\phi _t} \exp \left(\frac{3}{8}\sigma _s^2\right) \left[1+\mathrm{erf}\left(\frac{\sigma _s^2-s_{\rm crit}}{\sqrt{2\sigma _s^2}}\right)\right], \end{equation} \tag{ 41 }$$

which is identical to Equation (8) in Hennebelle & Chabrier (2011) and identical to the HC model, Equation (40), except that the critical density is defined according to either the KM, PN, or HC models. Thus, the multi-ff KM model is defined by using the threshold density s_crit = s_{crit, KM} from Equation (20) in the generalized solution of the multi-freefall SFR_ff, Equation (41).

2.4.5. The Multi-freefall PN Model

2.4.6. The Multi-freefall HC Model

After the detailed derivation of the six different SFR_ff models, we can now start to compare them. Figure 1 shows all six SFR_ff models: KM, PN, HC (left panels) and multi-ff KM, PN, HC (right panels) for a turbulent forcing parameter b = 0.4, corresponding to a statistical mixture of solenoidal and compressive modes in the turbulent forcing (Federrath et al. 2010b, Figure 8). When looking at the derivations above, it becomes clear that SFR_ff is a function of α_vir, $\mathcal {M}$, b, and β. The dependences enter through the definition of the critical densities in the different models, and through the variance of turbulent density fluctuations, Equation (4). We plot the analytically derived SFR_ff as a function of the virial parameter α_vir and the sonic Mach number $\mathcal {M}$ in each panel. Note that all these models are plotted for β → ∞, i.e., without taking magnetic fields into account yet. As shown in Table 1, each model has two fudge factors of order unity. The first one is 1/ϕ_t for all models (where the local efficiency was set to = 1 for simplicity in all models to facilitate the comparison), while the second one is ϕ_x, θ, and y_cut for the (multi-freefall) KM, PN, and HC models, respectively. We plot all curves for /ϕ_t = 1 to enable direct comparisons and used the favored values of the fudge factors by the different authors, ϕ_x = 1.12 (Krumholz & McKee 2005), θ = 0.35 (Padoan & Nordlund 2011), and y_cut = 0.1 (Hennebelle & Chabrier 2011).

Zoom In Zoom Out Reset image size

Dependence on α_vir. Let us first concentrate on the dependence of SFR_ff on the virial parameter. Since the virial parameter, Equations (15) and (16), is defined here as the ratio of twice the kinetic to the gravitational energy, it essentially measures how strongly the system is bound, and whether it is contracting (α_vir ≲ 1) or expanding (α_vir ≳ 1). Thus, we generally expect that the SFR should decrease with increasing α_vir because the cloud then becomes less bound and less likely to form stars. Indeed, the analytic SFR generally decreases with increasing α_vir in all models with the exception of the original PN model for which SFR_ff first increases for α_vir ≲ 1 and then decreases for α_vir ≳ 1. The increase comes from the freefall-time factor t_ff(ρ₀)/t_ff(ρ_crit) in the PN model, which leads to the factor exp (s_{crit, PN}/2) in Equation (32), and with the critical density from Equation (30) to SFR_ff∝α_vir for small α_vir. As expected though, this direct proportionality disappears in the multi-freefall PN model, as in the other two multi-freefall models (multi-ff KM and multi-ff HC).

Dependence on $\mathcal {M}$. The expected dependence of SFR_ff on the sonic Mach number is that SFR_ff should increase with increasing $\mathcal {M}$ because higher Mach number means stronger and denser local compression, leading to higher SFRs. Indeed, the Mach number dependence is generally similar in all models, i.e., SFR_ff increases with $\mathcal {M}$, with the exception of the original KM model, which has the weakest dependence on $\mathcal {M}$. For large α_vir, SFR_{ff, KM} increases, but only slowly, while for small α_vir, it stays constant or even decreases with increasing $\mathcal {M}$. Both the HC and multi-freefall HC models have the strongest positive correlation with the Mach number, such that for $\mathcal {M}\gtrsim 10$, SFR_{ff, HC} hardly depends on α_vir anymore (see also, Hennebelle & Chabrier 2011, Figure 1).⁷ The strong increase of SFR_ff with $\mathcal {M}$ in the two HC cases comes from the Mach number dependence of the thermal contribution to s_{crit, HC}, which is $\tilde{\rho }_{\rm crit,th} \propto \mathcal {M}^{-2}$, leading to a decreasing threshold density in the HC models, and thus to a higher SFR. This is the opposite compared to the KM and PN models, for both of which the critical density increases with the square of the Mach number (see Equations (20) and (30), respectively, or Table 1).

We also note the local minima of SFR_ff around $\mathcal {M}\approx 2$ in all models, except the HC and multi-freefall HC models. Those minima are spurious because they occur close to $\mathcal {M}=1$ for which the basic approach of shock-induced star formation must eventually break down as the system becomes transonic. Shocks require $\mathcal {M}>1$, by definition, but for rms Mach 1–2, a significant fraction of the system is transonic to subsonic. We thus conclude that all six models break down for the low Mach number regime, $\mathcal {M}\lesssim 2$. The rms sonic Mach numbers in real molecular clouds usually exceed unity by far (Larson 1981; Falgarone et al. 1992; Roman-Duval et al. 2011; Schneider et al. 2012), such that our analytic models are generally applicable to typical molecular clouds with $\mathcal {M}>2$.

Dependence on b. While the dependence on Mach number enters SFR_ff both through s_crit and σ²_s, the forcing dependence only enters through the forcing parameter b in σ²_s, Equation (4). Figure 2 shows SFR_ff as a function of the forcing parameter b for all models and three different Mach numbers ($\mathcal {M}=5$, 10, and 20). All curves are plotted for α_vir = 1, β → ∞, /ϕ_t = 1, and the standard fudge factors ϕ_x = 1.12, θ = 0.35, and y_cut = 0.1, respectively. We see that SFR_ff increases monotonically with b, from b = 1/3 (solenoidal forcing), over b = 0.4 (mixed forcing), to b = 1 (compressive forcing) in all models. This is expected because the density variance becomes larger for more compressive forcing, pushing a significant fraction of the gas to higher densities (Federrath et al. 2008b, 2010b; Konstandin et al. 2012a). Similar to the behavior with increasing Mach number, increasing the amount of direct compression induced by the turbulent forcing leads to higher local densities, and thus to higher SFRs with a typical increase of about an order of magnitude for compressive forcing compared to solenoidal forcing.

Zoom In Zoom Out Reset image size

Dependence on β. We expect that by adding magnetic energy to the system, the SFR should decrease because magnetic energy adds a stabilizing pressure to the system, counteracting gravitational collapse. Figure 3 shows the dependence of SFR_ff on plasma β in the six analytic models. We emphasize that only the original PN model had a magnetic-field dependence, coming from the dependence of s_{crit, PN} on β in Equation (30), and from the dependence of σ_s on β in Equation (4). However, we have extended all other analytic models (KM, HC, and multi-ff KM, PN, HC) to MHD, simply by applying the MHD version of σ_s, Equation (4) in all models, and replacing the sonic Mach number in the expressions for the critical density by the magnetic version $\mathcal {M}\rightarrow \mathcal {M}/\sqrt{1+1/\beta }$, introduced in Equation (19).

Zoom In Zoom Out Reset image size

As found in a detailed comparison of the analytically derived σ_s with numerical simulations of MHD turbulence in Molina et al. (2012), the standard-deviation–Mach number relation, Equation (4), breaks down for $\mathcal {M}_{\rm A}\lesssim 2$ because strongly sub-Alfvénic flows become highly anisotropic (e.g., Mac Low 1999; Cho & Vishniac 2000; Cho & Lazarian 2003; Beresnyak et al. 2005; Brunt et al. 2010; Esquivel & Lazarian 2011). Since the magnetic-field dependence of SFR_ff was introduced as an isotropic magnetic-pressure extension, the behavior of the analytic models for $\mathcal {M}_{\rm A}\lesssim 2$ is likely invalid. Thus, we only consider the trans- to super-Alfvénic regime with $\mathcal {M}_{\rm A}\gtrsim 2$. In this regime, SFR_ff decreases with increasing magnetic energy, i.e., decreasing β or $\mathcal {M}_{\rm A}$ in all models, as expected when adding a stabilizing magnetic pressure.

We use the adaptive mesh refinement (AMR; Berger & Colella 1989) code FLASH⁸ (Fryxell et al. 2000; Dubey et al. 2008) in version 2.5 to integrate the ideal, three-dimensional, MHD equations, including self-gravity,

$$\begin{eqnarray} & & \frac{\partial \rho }{\partial t} + \nabla \cdot \left(\rho \mathbf {v}\right)=0, \nonumber \\ & & \rho \left(\frac{\partial }{\partial t} + \mathbf {v}\cdot \nabla \right)\mathbf {v}= \frac{(\mathbf {B}\cdot \nabla)\mathbf {B}}{4\pi } - \nabla P_\star + \rho \left({\bf g} + {\bf F_{\rm stir}}\right), \nonumber \\ & & \frac{\partial E}{\partial t} + \nabla \cdot \left[\left(E+P_\star \right)\mathbf {v}- \frac{\left(\mathbf {B}\cdot \mathbf {v}\right)\mathbf {B}}{4\pi }\right] = \rho \mathbf {v}\cdot \left({\bf g}+{\bf F_{\rm stir}}\right), \nonumber \\ & & \frac{\partial \mathbf {B}}{\partial t} = \nabla \times \left(\mathbf {v}\times \mathbf {B}\right),\quad \nabla \cdot \mathbf {B}= 0, \end{eqnarray} \tag{ 42 }$$

where the gravitational acceleration of the gas g is the sum of the self-gravity of the gas and the contribution of sink particles (a subgrid model for collapse and accretion of star-forming regions in the simulations, explained in Section 3.3 below):

$$\begin{eqnarray} & & {\bf g} = -\nabla \Phi _{\rm gas} + {\bf g}_{\rm sinks}, \nonumber \\ & & \nabla ^2\Phi _{\rm gas} = 4\pi G\rho . \end{eqnarray} \tag{ 43 }$$

In the ideal MHD Equations (42), ρ, $\mathbf {v}$, $P_\star =P_{\rm th}+ 1/(8\pi)\left|\mathbf {B}\right|^2$, $\mathbf {B}$, and $E=\rho \epsilon _{\rm int} + (\rho /2)\left|\mathbf {v}\right|^2 + 1/(8\pi)\left|\mathbf {B}\right|^2$ denote gas density, velocity, pressure (thermal plus magnetic), magnetic field, and total energy density (internal plus kinetic, plus magnetic), respectively. The MHD equations are closed with a polytropic equation of state, $P_{\rm th}=c_{\rm s}^2\rho ^\Gamma$ with Γ = 1, such that the gas remains isothermal with a constant sound speed c_s = 0.2 km s⁻¹, corresponding to a temperature of T ≈ 11 K for gas with a mean molecular weight of 2.3. This is a reasonable approximation for dense, molecular gas of solar metallicity, over a wide range of densities (Wolfire et al. 1995; Omukai et al. 2005; Pavlovski et al. 2006; Glover & Mac Low 2007a, 2007b; Glover et al. 2010; Hill et al. 2011; Hennemann et al. 2012). Moreover, Glover & Clark (2012) find that the SFR is almost insensitive to the metallicity. Reducing the metallicity of the gas by two orders of magnitude reduces the time-averaged SFR by less than a factor of two. Thus, our conclusions remain intact, even though we neglect the detailed chemistry, cooling and heating processes in molecular clouds in this study.

We solve the MHD Equations (42) on three-dimensional, periodic grids with maximum resolutions of N³_res = 128³–1024³ grid points. These are all uniform-grid simulations, except for the N_res = 1024 simulation, where we use a root grid with 512³ cells and one level of AMR with a refinement criterion to ensure that the local Jeans lengths is covered with at least 32 grid cells, in order to resolve turbulent vorticity and magnetic-field amplification on the Jeans scale (Sur et al. 2010; Federrath et al. 2011c; Turk et al. 2012). We use a positive-definite MHD Riemann solver (Bouchut et al. 2007, 2010; Waagan 2009), which has been tested for efficiency, robustness, and accuracy in Waagan et al. (2011). This study shows that the MHD scheme keeps $\nabla \cdot \mathbf {B}$ errors at a negligible level and allows us to model extremely high Mach turbulence without producing unphysical states. This is particularly important for this study because we model supersonic turbulence on the largest scales of molecular clouds with rms Mach numbers as high as $\mathcal {M}\approx 50$ and compressive forcing, which produces density contrasts by several orders of magnitude, sometimes between two adjacent grid cells, because of multiple interactions of shocks and strong rarefaction waves, even before gravitational collapse sets in. Grid-based HD solvers often produce negative densities in such situations because of numerical post-shock oscillations. Such unphysical states are avoided by construction in the HLL3R Riemann scheme (Waagan et al. 2011) used here. The self-gravity of the gas, i.e., the gas–gas gravitational interaction (Equation (43)) is computed using a multi-grid Poisson solver (the FLASH2.5 version discussed in Ricker 2008), while the sink particle interactions are computed by direct N-body summation, as explained in Section 3.3 below. We note that the gravitational potential Φ_gas is computed with respect to the periodic boundary conditions specified in the simulations.

The ideal MHD Equations (42) do not contain any explicit kinematic viscosity and magnetic resistivity terms. However, any numerical scheme has an effective numerical viscosity ν and magnetic resistivity η due to the necessary discretization of the MHD equations. Even though the numerical viscosity depends on the specifications of the algorithm, it can be used to mimic the effects of explicit viscosity and resistivity (Benzi et al. 2008). It is important to realize, though, that the kinematic and magnetic Reynolds numbers that can be achieved with ideal MHD depend on the grid resolution. As shown in Federrath et al. (2011b), compressible, ideal MHD turbulence resolved with 128³ grid cells reaches kinematic Reynolds numbers Re = Lσ_V/ν ≈ 1500 and magnetic Reynolds numbers Rm = Lσ_V/η ≈ 3000. For Burgers (1948) scaling of the turbulence σ_v(ℓ)∝ℓ^1/2 (Equation (12) with p = 1/2), the Reynolds numbers scale as N^3/2_res as opposed to Kolmogorov (1941) scaling of the turbulence, σ_v(ℓ)∝ℓ^1/3 (Equation (12) with p = 1/3), leading to a Reynolds-number scaling of N^4/3_res. Thus, even in our highest resolution simulation with N_res = 1024, we only achieve Reynolds numbers Re ≈ (2.4–3.4) × 10⁴ and Rm ≈ (4.8–6.8) × 10⁴, depending on the scaling of the turbulence. In summary, although the flows we model exhibit fully developed turbulence (Frisch 1995), their Reynolds number are still considerably smaller than the ones inferred for real molecular clouds (see, e.g., Schober et al. 2012). We will thus study the resolution dependence of our results for the SFR below.

Previous numerical studies of non-driven turbulence have shown that supersonic turbulence decays in about a crossing time, irrespective of whether or not magnetic fields are included (Scalo & Pumphrey 1982; Mac Low et al. 1998; Stone et al. 1998; Mac Low 1999). The observed presence of turbulence has thus led to the conclusion that interstellar turbulence should be driven by some physical stirring mechanisms. Those mechanisms include supernova explosions and expanding, ionizing shells from previous cycles of star formation (McKee 1989; Krumholz et al. 2006; Balsara et al. 2004; Breitschwerdt et al. 2009; Peters et al. 2011; Goldbaum et al. 2011; Lee et al. 2012), gravitational collapse and accretion of material (Vazquez-Semadeni et al. 1998; Klessen & Hennebelle 2010; Elmegreen & Burkert 2010; Vázquez-Semadeni et al. 2010; Federrath et al. 2011c; Robertson & Goldreich 2012), and galactic spiral-arm compression of HI clouds (Dobbs & Bonnell 2008; Dobbs et al. 2008) and magnetorotational instability (MRI; Piontek & Ostriker 2007; Tamburro et al. 2009). On smaller scales, jets and outflows from young stellar objects have been suggested to drive turbulence (Norman & Silk 1980; Banerjee et al. 2007; Nakamura & Li 2008; Cunningham et al. 2009; Carroll et al. 2010; Wang et al. 2010). Turbulence in high-redshift galaxies is also likely driven by feedback from previous cycles of star formation (Green et al. 2010). A summary and comparison of driving mechanisms for interstellar turbulence is provided in Mac Low & Klessen (2004) and Elmegreen (2009). Mac Low & Klessen (2004) conclude that expanding shells are likely the dominant driver of interstellar turbulence in the star-forming parts of the Galaxy. More recently, Lee et al. (2012) also noted that the kinetic energy injected per unit of time by star-forming complexes via expansion of bubbles is about two-thirds of the luminosity required to maintain the observed velocity dispersions, supporting the view that expanding bubbles driven by massive star clusters from previous star formation are a major driver of turbulence in the Milky Way (see, e.g., the Cygnus X giant molecular cloud studied in Schneider et al. 2011).

It is important to realize that all these potential drivers (maybe with the exception of the MRI) are expected to primarily drive compressible modes in the velocity field, but do not directly excite solenoidal modes. However, even though the turbulence in molecular clouds might be driven compressively, solenoidal modes are indirectly excited by nonlinear interactions of multiple colliding shock fronts (Vishniac 1994; Sun & Takayama 2003; Kritsuk et al. 2007; Federrath et al. 2010b), by baroclinity, rotation and shear (Del Sordo & Brandenburg 2011), and by viscosity (Mee & Brandenburg 2006; Federrath et al. 2011b), such that supersonic turbulence driven by even purely compressive forcing contains about half of its kinetic power in solenoidal modes and the other half in compressible modes in the inertial range (Federrath et al. 2010b, Figure 14).

Modeling physical turbulent stirring mechanisms in numerical simulations requires assumptions about the spatial and temporal correlation of the turbulent forcing events. It is also still a matter of debate which of the physical mechanisms dominates the injection of turbulent energy on different cloud scales. Given these uncertainties, instead of trying to mimic one or more of the potential physical drivers of turbulence, we here use "driven turbulence in a box". From these simplified and idealized simulations, we can draw statistical conclusions about the role of turbulence for star formation, given the average simulation properties of a cloud (α_vir, $\mathcal {M}$, b, and β). In particular, our turbulent forcing approach allows us to evaluate the role of the mixture of velocity modes excited by a physical driver.

In practice, the stochastic forcing term F_stir is applied as a source term in Equations (42) to drive turbulence in the simulations. F_stir only contains large-scale modes, 1 < k < 3, where most of the power is injected at the k = 2 mode in Fourier space, which corresponds to half of the box size L in physical space. We thus model turbulent forcing on large scales, as favored by molecular cloud observations (e.g., Ossenkopf & Mac Low 2002; Heyer et al. 2006; Brunt et al. 2009; Gaensler et al. 2011; Roman-Duval et al. 2011). Smaller scales, k > 3, are not affected directly by the forcing, such that turbulence can develop self-consistently on these scales. We use the Ornstein–Uhlenbeck (OU) process to model F_stir, which is a well-defined stochastic process with a finite autocorrelation timescale (Eswaran & Pope 1988; Schmidt et al. 2006), leading to a smoothly varying stochastic force field in space and time. Details about the OU process and the forcing applied in this study can be found in Schmidt et al. (2009), Federrath et al. (2010b), and Konstandin et al. (2012a). However, the essential point of our forcing approach is that we can adjust the mixture of solenoidal and compressive modes of F_stir. This is achieved by decomposing a given vector field with random mixtures into its solenoidal and compressive parts, by applying the projection tensor $\mathcal {\underline{P}}^{\,\zeta }({\mathbf {k}})$ in Fourier space. In index notation, this tensor reads

$$\begin{equation} \mathcal {P}_{ij}^\zeta = \zeta \,\mathcal {P}_{ij}^\perp +(1-\zeta)\,\mathcal {P}_{ij}^\parallel = \zeta \,\delta _{ij}+(1-2\zeta)\,\frac{k_i k_j}{|k|^2}, \end{equation} \tag{ 44 }$$

where δ_ij is the Kronecker symbol, and $\mathcal {P}_{ij}^\perp = \delta _{ij} - k_i k_j / k^2$ and $\mathcal {P}_{ij}^\parallel = k_i k_j / k^2$ are the solenoidal and compressive projection operators, respectively. The ratio of compressive power to total power in F_stir can be derived from Equation (44) by evaluating the norm of the compressive component of the projection tensor and dividing it by the total injected power, resulting in

$$\begin{equation} \frac{F_{\rm comp}}{F_{\rm tot}} = \frac{(1-\zeta)^2}{1-2\zeta +3\zeta ^2}\,, \end{equation} \tag{ 45 }$$

for three-dimensional space (Schmidt et al. 2009; Federrath et al. 2010b). The projection operator serves to construct a purely solenoidal force field by setting ζ = 1, while for ζ = 0 a purely compressive force field is obtained. Any combination of solenoidal and compressive modes can be constructed by choosing ζ ∈ [0, 1]. Here we compare simulations with ζ = 1 (sol.), ζ = 1/2 (mix.), and ζ = 0 (comp.). A detailed study of the forcing dependence of the b-parameter entering the expression for the variance of the density PDF, Equations (4) and (5), is provided in Federrath et al. (2010b, Figure 8), where they measure b as a function of the forcing parameter ζ.

In order to model collapse and accretion of star-forming gas in the simulations, we use a subgrid model called "sink particles," which is a method originally invented by Bate et al. (1995) for smoothed particle hydrodynamics, and first adopted for Eulerian AMR simulations by Krumholz et al. (2004). In Krumholz et al. (2004), a Lagrangian sink particle is introduced, if the gas reaches a given density. However, sink particles are supposed to represent bound objects that are going into collapse, and thus, a density threshold as the only criterion for sink particle creation is insufficient (Federrath et al. 2010a). Based on the ideas of Bate et al. (1995) and Krumholz et al. (2004), we use an advanced AMR-based approach for sink particles in which only bound and collapsing gas is accreted, thus avoiding the creation of spurious sink particles (for a detailed analysis, see Federrath et al. 2010a). The key feature of this approach is to define a control volume around cells that exceed the density threshold set by the resolution criterion to avoid artificial fragmentation. Truelove et al. (1997) found that the Jeans length must be resolved with at least four grid cells to avoid artificial fragmentation, leading to a resolution-dependent density threshold criterion for the creation of sink particles:

$$\begin{equation} \rho _{\rm sink}= \frac{\pi c_{\rm s}^2}{4G\,r_{\rm sink}^2}\,, \end{equation} \tag{ 46 }$$

where the sink particle accretion radius r_sink is set to 2.5 grid-cell lengths at the maximum level of refinement, corresponding to half a Jeans length at ρ_sink, such that the Jeans length is still resolved with five grid cells prior to potential sink particle creation to avoid artificial fragmentation. Grid cells exceeding the density threshold given by Equation (46), however, do not form sink particles right away. First, a spherical control volume with radius r_sink is defined around the cell exceeding ρ_sink within which additional checks for gravitational instability and collapse are performed. We check whether the gas

• 1.  
  is on the highest level of refinement,
• 2.  
  is converging from all directions in the rest frame of the central cell (negative radial velocity),
• 3.  
  is at a local gravitational potential minimum,
• 4.  
  is bound (|E_grav| > E_kin + E_th + E_mag),
• 5.  
  is Jeans unstable, and
• 6.  
  is not within r_sink of an existing sink particle.

If all these checks are passed, a sink particle is created in the center of the control volume (see Federrath et al. 2010a). This procedure avoids spurious sink particle formation and allows us to trace only truly collapsing and star-forming gas. Given the checks above, it is clear that in some cases, a sink particle is not necessarily formed even though the density threshold is exceeded. This does not mean, however, that such gas would be subject to artificial gravitational fragmentation. Since the checks did not allow sink particle creation, the gas in the control volume was not collapsing and/or not bound, there is no need to worry about artificial fragmentation at this stage, even though the density threshold was exceeded. This can happen quite frequently in supersonic turbulence because shocks can push the gas density above the threshold, even though this gas is not necessarily gravitationally bound after the shock passage.

Once a sink particle is created, it can gain mass by accreting gas from the AMR grid, but only if this gas exceeds the threshold density, is inside the sink particle accretion radius, is bound to the particle, and is collapsing toward it. If all these criteria are fulfilled, the excess mass above the density threshold defined by Equation (46) is removed from the MHD system and added to the sink particle, such that mass, momentum and angular momentum are conserved by construction (see Federrath et al. 2010a, 2011a, for details).

All contributions to the gravitational interactions between the gas on the grid and the sink particles are computed by direct N-body summation over all grid cells and sink particles (gas–sink, sink–gas, and sink–sink), using gravitational spline softening inside the sink particle radius to avoid singularities during close encounters. The softening only affects scales that are below the grid-resolution cutoff set by the sink particle accretion radius. A second-order accurate Leapfrog integrator is used to advance the sink particles on a time step that allows us to resolve close and highly eccentric orbits of sink particles without introducing significant errors on super-resolution grid scales.

Starting from a uniform density distribution and zero velocities, the forcing term F_stir in Equations (42) excites turbulent motions. First, we evolve the MHD equations for two turbulent crossing times, $2T=L/(\mathcal {M}c_{\rm s})$ without self-gravity, in order to establish fully developed, compressible turbulence (e.g., Klessen et al. 2000; Klessen 2001; Heitsch et al. 2001; Federrath et al. 2009, 2010b; Price & Federrath 2010; Micic et al. 2012). We do not include the gravity terms until t = 2 T, in order to avoid our measurements of the SFR being contaminated by this rather artificial initial transient phase, during which the system is building up a turbulent cascade (Schmidt et al. 2009). After that we solve the full system of MHD equations (42) and (43) including self-gravity and formation of sink particles. For practical purposes, we reset the time t = 2 T to t = 0 t_ff(ρ₀), which is the time when turbulence is fully established and star formation is allowed to proceed. We note that this procedure is slightly different from setting up a simulation with power-law velocity scaling drawn from Gaussian random seeds as an initial condition, commonly applied in numerical star formation studies (e.g., Bate et al. 2003; Clark et al. 2005; Krumholz et al. 2007; Price & Bate 2008, 2009; Smith et al. 2008; Federrath et al. 2010a; Walch et al. 2010; Girichidis et al. 2011). In those cases, the initial random velocity field is imposed on top of a given density profile (often constant density or radial power-law distributions), such that density and velocity fields have no causal connection. Here, the initial density and velocity fields at t = 0 are consistently coupled via the equations of (magneto)hydrodynamics. We also keep driving the turbulence instead of imposing only an initial Gaussian perturbation as in the studies mentioned above.

All our numerical simulations and their basic parameters are listed in Table 2. Each model has a unique name, starting with "GT" (for "GravTurb"), followed by the maximum grid resolution ("128," "256," "512," and "1024"), the forcing type ("s": solenoidal; "m": mixed; and "c": compressive), and the Mach number ("M3," "M5," "M10," "M20," and "M50"). Models with an initially uniform magnetic field in the z-direction through the simulation box are additionally denoted "B1," "B3," and "B10," corresponding to B₀ = 1, 3, and 10 μG, respectively. Different random sequences with the same statistical properties for the turbulent forcing are indicated by "(s1)," "(s2)," and "(s3)" at the end of the model name, indicating that random "(seed1)," "(seed2)," or "(seed3)" was used. Columns 2–10 in Table 2 list the maximum numerical resolution, type of forcing, mean density ρ₀, box size L, the total mass M_c, large-scale velocity dispersion σ_V, initial magnetic-field strength B₀, initial plasma β₀, and the virial parameter α_vir, ° computed with Equation (15).

Table 2. Basic Parameters of the Numerical Models of Forced, Supersonic, Self-gravitating, (Magneto)hydrodynamic Turbulence

Model	Nres	Forcing	ρ0	L	Mc	σV	B0	β0	αvir, °	αvir	$\mathcal {M}$	b	β	$\mathcal {M}_{\rm A}$
			(g cm−3)	(pc)	(M☉)	(km s−1)	(μG)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)	(14)	(15)
(01) GT256sM3	256	sol	5.8 × 10−19	3.3 × 10−1	3.1 × 102	0.59	0	∞	0.07	1.4	2.9	1/3	∞	∞
(02) GT512sM3	512	sol	5.8 × 10−19	3.3 × 10−1	3.1 × 102	0.59	0	∞	0.07	1.4	3.0	1/3	∞	∞
(03) GT256mM3	256	mix	5.8 × 10−19	3.3 × 10−1	3.1 × 102	0.61	0	∞	0.08	1.1	3.1	0.4	∞	∞
(04) GT256cM3	256	comp	5.8 × 10−19	3.3 × 10−1	3.1 × 102	0.58	0	∞	0.07	0.46	2.9	1	∞	∞
(05) GT512cM3	512	comp	5.8 × 10−19	3.3 × 10−1	3.1 × 102	0.58	0	∞	0.07	0.48	2.9	1	∞	∞
(06) GT256sM5	256	sol	3.3 × 10−21	2.0 × 100	3.9 × 102	1.0	0	∞	1.0	8.0	5.0	1/3	∞	∞
(07) GT256mM5	256	mix	3.3 × 10−21	2.0 × 100	3.9 × 102	0.99	0	∞	0.98	5.4	5.0	0.4	∞	∞
(08) GT256cM5	256	comp	3.3 × 10−21	2.0 × 100	3.9 × 102	0.91	0	∞	0.82	1.5	4.5	1	∞	∞
(09) GT128sM10	128	sol	8.2 × 10−22	8.0 × 100	6.2 × 103	2.1	0	∞	1.1	11.	10.	1/3	∞	∞
(10) GT256sM10	256	sol	8.2 × 10−22	8.0 × 100	6.2 × 103	2.1	0	∞	1.1	12.	10.	1/3	∞	∞
(11) GT512sM10	512	sol	8.2 × 10−22	8.0 × 100	6.2 × 103	2.1	0	∞	1.1	12.	10.	1/3	∞	∞
(12) GT512mM10 (s1)	512	mix	8.2 × 10−22	8.0 × 100	6.2 × 103	2.1	0	∞	1.1	4.5	11.	0.4	∞	∞
(13) GT512mM10B1 (s1)	512	mix	8.2 × 10−22	8.0 × 100	6.2 × 103	2.1	1	8.2	1.1	5.4	10.	0.4	2.8	12.
(14) GT512mM10 (s2)	512	mix	8.2 × 10−22	8.0 × 100	6.2 × 103	2.2	0	∞	1.2	8.4	11.	0.4	∞	∞
(15) GT512mM10B1 (s2)	512	mix	8.2 × 10−22	8.0 × 100	6.2 × 103	2.2	1	8.2	1.2	9.5	11.	0.4	1.8	10.
(16) GT256mM10 (s3)	256	mix	8.2 × 10−22	8.0 × 100	6.2 × 103	2.0	0	∞	1.0	5.9	10.	0.4	∞	∞
(17) GT512mM10 (s3)	512	mix	8.2 × 10−22	8.0 × 100	6.2 × 103	2.0	0	∞	1.0	5.9	10.	0.4	∞	∞
(18) GT512mM10B1 (s3)	512	mix	8.2 × 10−22	8.0 × 100	6.2 × 103	2.0	1	8.2	0.97	6.4	9.9	0.4	3.6	13.
(19) GT256mM10B3 (s3)	256	mix	8.2 × 10−22	8.0 × 100	6.2 × 103	1.8	3	0.92	0.81	8.4	9.0	0.4	0.20	2.9
(20) GT512mM10B3 (s3)	512	mix	8.2 × 10−22	8.0 × 100	6.2 × 103	1.8	3	0.92	0.83	8.7	9.1	0.4	0.18	2.7
(21) GT256mM10B10 (s3)	256	mix	8.2 × 10−22	8.0 × 100	6.2 × 103	1.8	10	0.08	0.79	6.6	8.9	0.4	0.04	1.3
(22) GT128cM10	128	comp	8.2 × 10−22	8.0 × 100	6.2 × 103	1.8	0	∞	0.81	1.2	9.0	1	∞	∞
(23) GT256cM10	256	comp	8.2 × 10−22	8.0 × 100	6.2 × 103	1.8	0	∞	0.85	1.1	9.2	1	∞	∞
(24) GT512cM10	512	comp	8.2 × 10−22	8.0 × 100	6.2 × 103	1.9	0	∞	0.87	1.1	9.4	1	∞	∞
(25) GT256sM20	256	sol	2.1 × 10−22	3.2 × 101	9.9 × 104	4.1	0	∞	1.0	11.	20.	1/3	∞	∞
(26) GT256mM20	256	mix	2.1 × 10−22	3.2 × 101	9.9 × 104	4.2	0	∞	1.1	4.5	21.	0.4	∞	∞
(27) GT256cM20	256	comp	2.1 × 10−22	3.2 × 101	9.9 × 104	4.0	0	∞	1.0	0.60	20.	1	∞	∞
(28) GT256sM50	256	sol	3.3 × 10−23	2.0 × 102	3.9 × 106	10.	0	∞	1.1	12.	52.	1/3	∞	∞
(29) GT512sM50	512	sol	3.3 × 10−23	2.0 × 102	3.9 × 106	10.	0	∞	1.1	13.	52.	1/3	∞	∞
(30) GT256mM50	256	mix	3.3 × 10−23	2.0 × 102	3.9 × 106	10.	0	∞	1.0	7.0	51.	0.4	∞	∞
(31) GT512mM50	512	mix	3.3 × 10−23	2.0 × 102	3.9 × 106	10.	0	∞	1.1	7.4	51.	0.4	∞	∞
(32) GT256cM50	256	comp	3.3 × 10−23	2.0 × 102	3.9 × 106	9.8	0	∞	0.95	0.54	49.	1	∞	∞
(33) GT512cM50	512	comp	3.3 × 10−23	2.0 × 102	3.9 × 106	9.9	0	∞	0.99	0.56	50.	1	∞	∞
(34) GT1024cM50	1024	comp	3.3 × 10−23	2.0 × 102	3.9 × 106	10.	0	∞	1.00	0.55	50.	1	∞	∞

Notes. Column 1: simulation name. Columns 2–10: maximum grid resolution in one direction of the three-dimensional cubic domain, mode of forcing (solenoidal, mixed, compressive), mean density, linear box size, total mass, velocity dispersion on the box scale, mean magnetic-field strength (in the z-direction of the domain), initial plasma β₀, and virial parameter based on Equation (15). Columns 11–15: time-averaged virial parameter based on Equation (16), computed directly from the three-dimensional gas distribution, the sonic Mach number, forcing parameter, ratio of thermal to magnetic pressure (plasma β), and Alfvén Mach number. To guide the eye, horizontal lines separate models with different sonic Mach numbers.

Download table as:  ASCIITypeset image

Columns 11–15 are derived quantities, measured as space and time averages after turbulence is fully established, t ⩾ 0, until 20% of the original cloud mass is accreted onto sink particles, i.e., the star formation efficiency has reached SFE = 20%. We list the average virial parameter α_vir, the sonic Mach number $\mathcal {M}$, forcing parameter b, plasma β, and Alfvén Mach number $\mathcal {M}_{\rm A}$. The instantaneous virial parameter, Equation (16), in Column 11 of Table 2 is computed as α_vir = 2E_kin/|E_grav| = ∑M_iv²_i/|∑M_iΦ_{gas, i}| from the gravitational potential Φ_gas returned by the Poisson solver (see Section 3.1), as a sum over all grid cells i with mass M_i and velocity v_i. We note that this is different from the value α_vir, ° obtained from Equation (15) and listed in Column 10, which assumes a homogenous, spherical density distribution. In contrast, we obtain highly inhomogeneous density distributions in our compressible, turbulent clouds. We thus prefer to compute α_vir based on the three-dimensional density field as explained above.⁹ By analogy, the sonic and Alfvén Mach numbers, as well as β are computed as spatial rms averages over all cells in the simulation box as a function of time, followed by averaging over time. We will show in the next section that this approach is justified because we find that all those parameters do not vary significantly with time during star formation. The value of the forcing parameter b was not determined by averaging because it was already measured in Federrath et al. (2010b, Figure 8), giving best-fit values b = 1/3, 0.4, and 1 for solenoidal, naturally mixed, and compressive forcing of the turbulence, respectively.

We do not include any data or discussion of the state of the clouds after SFE = 20% is reached because at that point in time local feedback processes would have likely altered the subsequent evolution of the clouds so drastically that we cannot trust our results for higher SFE. Even before that, inclusion of feedback processes might change the results, at least locally. For example, we expect the amount of accreted gas to be reduced, if feedback were included (e.g., Wang et al. 2010; Peters et al. 2011). This fact can be accounted for by adjusting the local efficiency parameter introduced in Equation (7) to values < 1 for all the models discussed here. We return to this issue when we compare our simulations with the observational data in Section 6.

The basic model parameters in Table 2 were chosen to roughly follow observed properties of molecular clouds, covering a range of cloud sizes L ≈ 0.3–200 pc, masses M_c ≈ 300 to 4 × 10⁶ M_☉, and velocity dispersions σ_V ≈ 0.6–10 km s⁻¹ (e.g., Larson 1981; Solomon et al. 1987; Falgarone et al. 1992), with typical cloud scalings summarized and discussed in Mac Low & Klessen (2004) and McKee & Ostriker (2007). However, even though most real clouds may roughly follow such an average scaling, the scatter around that average is typically about an order of magnitude or more in terms of mass, density, and velocity dispersion for a given cloud size (e.g., Heyer et al. 2009; Roman-Duval et al. 2010). The procedure used here to determine the initial cloud parameters in the simulations is as follows. First, for a given target Mach number, we determine the appropriate size of the cloud by inverting the observed velocity dispersion–size relation given by Equation (12). Having the size and velocity dispersion, we then set the virial parameter given by Equation (15) to a value close to unity. The only exceptions are the $\mathcal {M}\sim 3$ models, where we set it to α_vir, ° ≈ 0.07, because this turned out to give actual virial parameters α_vir closer to unity after the turbulence had been fully established (compare Columns 10 and 11 in Table 2). Using the initial guess of α_vir, °, we then solve for the mass of the cloud by inverting Equation (15). From the mass and size, we compute the mean density of the model cloud.

It is important to note that the actual virial parameter obtained after two turbulent crossing times can be up to an order of magnitude different from the initial guess provided by Equation (15), depending on the Mach number and forcing of the model (see Table 2). This is because the density distribution in the state of fully developed supersonic turbulence is highly inhomogeneous and is not well described by Equation (15). Thus, we do not know the virial parameter that arises in the regime of fully developed turbulence a priori. The α_vir in the turbulent phase is typically higher (except for the compressive-forcing cases at high Mach numbers, $\mathcal {M}\sim 20$ and 50) than the one computed in Equation (15), because we also use periodic boundary conditions. This reduces the gravitational binding energy of the system compared to an isolated system (as assumed in Equation (15)). Real clouds are neither periodic nor isolated, but using periodic boundaries, we mimic the effects of the surrounding medium on the region studied in our computational boxes (discussed further in Section 7). We emphasize that the virial parameters obtained here are consistent with observations, given that observational estimates of α_vir are usually obtained based on Equation (15) or column-density versions of it.

Magnetic-field strengths for the MHD simulations were chosen to be consistent with the range observed in clouds (e.g., Crutcher 1999; Crutcher et al. 2010). We vary the magnetic field for simulations with mixed forcing and fixed sonic Mach number $\mathcal {M}\sim 10$, which gives us a good indication of the role of magnetic fields for typical molecular cloud properties. Heiles & Troland (2005) and Crutcher et al. (2010) show that most clouds with number densities in the range 10–10⁴ cm⁻³ have magnetic-field strengths in the range B_z ≈ 1–10 μG, with an apparent peak of the distribution at around B_z ≈ 3 μG. Our MHD simulations have mean densities of about 200 cm⁻³, so we decided to compare models with B_z = 1, 3, and 10 μG, in order to cover the observed range of line-of-sight magnetic-field strengths.

4.1.1. Effects of the Magnetic Field

Figure 4 shows the time evolution of column density snapshots (from top to bottom) for models with mixed forcing at $\mathcal {M}=10$ and 512³ resolution for initial magnetic fields B₀ = 0, 1, and 3 μG (left, middle, and right panels). Key initial parameters (box size, total mass, etc.), the time in units of t_ff(ρ₀), the SFE, and the number of sink particles formed are given in each panel. The top row shows the gas at t = 0, i.e., when turbulence is fully developed and self-gravity is switched on. We see shocks and large-scale structure induced by the large-scale turbulence with column density contrasts ranging over more than four orders of magnitude. Comparing the purely HD run (left) with the two magnetized runs (middle and right), we see that shocks become smoother and density contrasts slightly decrease as the magnetic-field strength increases. This is because magnetic fields act like a cushion, reducing density fluctuations, due to the additional magnetic pressure parameterized either by plasma β or the Alfvén Mach number $\mathcal {M}_{\rm A}$ (see Equation (4) or (5) and Molina et al. 2012), the time-averaged values of which are given in Table 2. At later times, the gas starts collapsing locally at sites previously compressed by the supersonic turbulence, at which point local filaments become more and more massive as they accrete gas from the surrounding and eventually become so dense that these cores have to be replaced with sink particles, allowing us to advance the simulations to later times (see Section 3.3). The radius r_sink of the sink particles is determined by the numerical resolution constraint and is given in each panel, as soon as sink particles have formed. Our resolution is insufficient to resolve individual stars, but the sink particles can be regarded as dense, bound cores in our simulations.

Zoom In Zoom Out Reset image size

Comparing the runs with different magnetic-field strengths in Figure 4, we see two important effects with increasing magnetic field: (1) a reduction of fragmentation, i.e., fewer sink particles have formed by the end of the simulations at SFE = 20% and (2) reaching a given SFE takes longer, i.e., the core formation rate and hence the SFR are reduced. For instance, when the SFE has reached 20%, the runs with B₀ = 0, 1, and 3 μG have formed 109, 71, and 63 sink particles in 0.71, 0.85, and 1.1 t_ff(ρ₀), respectively.

The higher the magnetic field, the larger the topologically connected structures, compared to the more fragmented and dispersed filaments in the purely hydrodynamical run. Comparing numerical simulations and observations with filament-tracking tools (e.g., André et al. 2010; Men'shchikov et al. 2010; Arzoumanian et al. 2011; Hill et al. 2011; Schneider et al. 2012) or polarization analyses (e.g., Burkhart et al. 2012) may eventually help to reveal the role of magnetic fields. In particular, the orientation of magnetic fields might tell us about its dynamical influence (Schneider et al. 2010; Li & Henning 2011; Peretto et al. 2012). In Figure 5, we show the time evolution of column density snapshots with local magnetic-field vectors computed by a mass-weighted average along the line of sight superimposed, for the run with B₀ = 3 μG for t/t_ff = 0 (top) and SFE = 10% (bottom). The magnetic field grows due to compression of the field lines and due to dynamo action (Sur et al. 2010; Federrath et al. 2011c; Bertram et al. 2012), particularly in regions where dense cores accumulate and form clusters. The magnetic field is very intermittent and shows no particularly preferred direction in the cluster centers because the gas motions are so chaotic that the magnetic-field direction changes frequently. The magnetic field is of moderate strength compared to the turbulence in this case, shown by the average super-Alfvénic Mach number in this simulation, $\mathcal {M}_{\rm A}\approx 2.7$ (see Table 2). The field strengths are consistent with observations in typical molecular clouds. On scales larger than molecular clouds and on Galactic scales though, the turbulence might be trans-Alfvénic rather than super-Alfvénic, which would naturally lead to more aligned magnetic-field structures there (e.g., Beck et al. 1996; Heiles & Troland 2005; Li & Henning 2011).

Zoom In Zoom Out Reset image size

4.1.2. Effects of Turbulent Forcing and Sonic Mach Number

After having looked at the time evolution of column density snapshots in mixed, $\mathcal {M}\sim 10$ simulations, we now focus on the gas morphology when SFE = 10%, representing a typical molecular cloud value, but comparing different forcing and sonic Mach numbers. Figure 6 shows column density projections of the 512³ runs with solenoidal forcing (left panels) and compressive forcing (right panels) at $\mathcal {M}\sim 3$ (top), $\mathcal {M}\sim 10$ (middle), and $\mathcal {M}\sim 50$ (bottom). Note the different length and mass scales probed in these images, with box sizes of L = 0.3, 8 and 200 pc, and masses of M_c = 310, 6.2 × 10³ and 3.9 × 10⁶ M_☉, respectively. Since the resolution is fixed, the sink particle radii vary from r_sink = 335 AU over r_sink = 0.04 pc, up to 1 pc. Thus, neither of those represents stars, but rather star clusters in the largest-scale runs and potentially protostellar accretion envelopes in the smallest-scale runs. The scale and mass sequence from the bottom to the top panels in Figure 6 can be interpreted as zooms into patches of larger-scale runs and re-simulating these patches with higher resolution in successively smaller boxes. Clearly, these images emphasize how artificial this kind of numerical experiment is, yet real molecular clouds exhibit similar hierarchical structures (Falgarone et al. 1992; Ossenkopf & Mac Low 2002), often characterized as fractals (Scalo 1990; Elmegreen & Falgarone 1996; Stutzki et al. 1998; Sánchez et al. 2005; Roman-Duval et al. 2010). The fractal dimension D inferred from different techniques (Δ-variance, box counting, mass–size relation, and perimeter-area method) was shown to vary between D ≈ 2.6 for purely solenoidal and D ≈ 2.3 for purely compressive forcing, in the range of observational determinations (Federrath et al. 2009), and is consistent with theoretical ideas to explain the slope of the stellar IMF (Chabrier & Hennebelle 2011). As can be seen in Figure 6, compressive forcing produces more sheet-like structures (planar shocks), while solenoidal forcing produces more volume-filling structures, providing a visual explanation for the dependence of D on the forcing.

Zoom In Zoom Out Reset image size

Besides the morphological distinctions, the most striking difference between solenoidal and compressive forcings is the timescale of core and star formation (compare t/t_ff in the upper right corner of each panel in Figure 6). For fixed Mach number, cloud size, and mass, compressive forcing accelerates the conversion of gas into stars compared to solenoidal forcing by factors of 4, 8, and 12 for the $\mathcal {M}\sim 3$, 10, and 50 runs, respectively, when SFE ∼ 10%. This result emphasizes the important role of the turbulent forcing for setting the SFR.

4.2.1. Effects of the Forcing, Random Seed, and Resolution

We now turn to the detailed time evolution and determination of the SFR in the simulations. Figure 7 shows the time evolution of dynamical quantities $\mathcal {M}$, α_vir, and SFE for models with $\mathcal {M}\sim 10$ and B₀ = 0 for solenoidal and compressive forcing at numerical resolutions of 128³, 256³, and 512³ grid cells, and for mixed forcing with three different random seeds of the turbulent forcing. The Mach number (top panel) shows variations of order 10% around the target Mach number of $\mathcal {M}\sim 10$ for each simulation, and some systematic variations with different forcing and random seeds. The differences between solenoidal, mixed, and compressive forcings are caused by stronger dissipation with more compressive forcing, requiring a higher forcing amplitude to reach the same Mach number than in solenoidal forcing. We adjust the amplitude of the forcing such that the gas reaches a given Mach number in the fully developed turbulent phase. Since the value of $\mathcal {M}$ depends on nonlinear dissipation properties, i.e., strengths of shocks and amount of vorticity generated, and thus on the Mach number of the turbulence (Federrath et al. 2011b), the time-averaged rms Mach number for a given forcing amplitude cannot be predicted a priori and must be adjusted iteratively by running test simulations with different forcing amplitude and measuring the time-averaged rms Mach number, resulting in some deviation of the actual Mach number from the target Mach number (see the time-averaged $\mathcal {M}$ for each model in Table 2). The temporal fluctuations and the differences between random seeds, however, are purely statistical. In order to compare our simulation data with the analytic theories, we thus always use the volume- and time-averaged quantities entering the theoretical models from Section 2.5.

Zoom In Zoom Out Reset image size

The middle panel of Figure 7 shows α_vir(t). As for $\mathcal {M}(t)$, the resolution dependence is only marginal and significantly less than the statistical fluctuations (see also Kitsionas et al. 2009; Price & Federrath 2010; Kritsuk et al. 2011a, showing that one-point statistics are typically well converged with grid resolutions of 256³ cells). This demonstrates that the length scales of the dominant gravitational structures are resolved well enough in our present numerical experiments and that our definition of α_vir is robust with respect to changes in the numerical resolution.

The difference of α_vir ≡ 2 E_kin/|E_grav| between the forcings deserves some attention. While the $\mathcal {M}\sim 10$ runs with solenoidal forcing have α_vir ≈ 12, the compressive ones have α_vir ≈ 1.1 (see Table 2), even though the Mach number is similar and the mass of the clouds is identical. In Figure 6, we saw that compressive forcing produces more locally compressed structures than solenoidal forcing, resulting in an overall higher gravitational binding energy |E_grav| compared to solenoidal forcing. The total kinetic energy E_kin on the other hand is the same within a factor of ∼2, which means that the factor of ∼10 difference in α_vir is primarily due to the difference in |E_grav|. This shows that comparing simple theoretical estimates of the virial parameter, solely based on the total mass as a measure for |E_grav| (as, e.g., assumed in Equation (15)), should be considered with great caution because such an estimate ignores the internal structure of the clouds. Thus, we prefer to estimate α_vir based on the actual spatial distribution, as we have done in Figure 7 and in Table 2 for all models. We emphasize that this direct comparison of α_vir, ° with α_vir performed here means that observational estimates of the virial parameter based on global measures such as described by Equation (15) or alike are only accurate within an order of magnitude. Measurements of gravitational (in)stability based on the actual column density distribution of filaments in Herschel observations of the Gould Belt (GB) by André et al. (2010), for example, are thus likely more accurate and meaningful than estimates based on uniform-density, spherical approximations such as Equation (15).

The bottom panel of Figure 7 shows the time evolution of the total mass accreted by sink particles, divided by the total cloud mass, i.e., the SFE. We measure the slope of these curves by fitting a linear function in the interval SFE = [4% , 20%], which gives the SFR_ff for each model quoted in the legend. We choose to set the lower limit of the fit range to SFE = 4% because the initial accretion phase is highly nonlinear with a fast increase in slope, after which the accretion becomes roughly linear in time, such that the slope is reasonably well defined for most of the models.

First, we study the dependence of SFR_ff on the random seed. The three models with mixed forcing and different random seeds (seed1, seed2, seed3) exhibit variations in SFR_ff by a factor of 1.5. However, other seeds might deviate further from this, such that the factor 1.5 in SFR_ff is a lower limit for the uncertainty introduced by the random seed. When we compare mixed-forcing models at $\mathcal {M}\sim 10$ with different magnetic-field strengths later, we always compare runs with seed3 because the SFR_ff for seed3 is in between the ones measured for seed1 and seed2, thus giving the best average behavior for the data at hand.

Finally, we investigate the resolution dependence of our simulations with solenoidal and compressive forcings in Figure 7. For resolutions of 128³, 256³, and 512³ grid cells, we find that SFR_ff = 0.27, 0.17, and 0.14 for solenoidal forcing, and SFR_ff = 1.58, 2.27, and 2.75 for compressive forcing, respectively. Thus, with increasing resolution, SFR_ff is decreasing for solenoidal forcing, but increasing for compressive forcing. The difference in SFR_ff between 128³ and 256³ is a factor of 0.63 for solenoidal forcing, and a factor of 1.44 for compressive forcing. These factors become smaller when we compare the 256³ with the 512³ simulations, giving factors of 0.82 for solenoidal forcing and 1.21 for compressive forcing. Thus our results converge with increasing resolution. Moreover, we can estimate SFR_ff in the limit of infinite resolution from extrapolating the convergence behavior. Doing this, we see that our measurements of SFR_ff at 128³ resolution are converged only within a factor of about 2.5, so we discard the two 128³ simulations (GT128sM10 and GT128cM10) in all of the following. In contrast, the 256³ data are converged to within a factor of 1.5 for both solenoidal and compressive forcings, which is similar to the uncertainty introduced by varying the random seed as discussed in the previous paragraph. Thus, differences larger than a factor of 1.5 in SFR_ff between models with different physical parameters are likely of physical rather than numerical or statistical origin. For instance, the SFR_ff for compressive forcing with $\mathcal {M}\sim 10$ is more than an order of magnitude larger than the SFR_ff of the respective solenoidal simulation, demonstrating the physical importance of the turbulent forcing for controlling the SFR.

4.2.2. Effects of Increasing the Sonic Mach Number

Having looked at models with different forcing, random seed, and resolution, we now study models with varying Mach number. Figure 8 shows the same as Figure 7, but for compressive-forcing models with Mach $\mathcal {M}\sim 3$, 5, 10, 20, and 50. The Mach number increases slightly with time in all models, which is caused by local accelerations during collapse. This only accounts for a few percent at most. The exception is the $\mathcal {M}\sim 3$ model for which $\mathcal {M}$ increases by almost a factor of two (top panel). This is because the $\mathcal {M}\sim 3$ model is more gravitationally unstable initially, indicated by the virial parameter (middle panel), which increases to around unity, similar to the other models. The SFR_ff generally increases with Mach number due to the stronger local compressions created at higher $\mathcal {M}$ (bottom panel). The only exception is again the $\mathcal {M}\sim 3$ model, which has a slightly higher SFR_ff than the $\mathcal {M}\sim 5$ model, because the $\mathcal {M}\sim 3$ model is more unstable and starts collapsing globally, while this is not the case in the other models. The difference of SFR_ff between $\mathcal {M}\sim 5$ and 50 is about a factor of 4.4.

Zoom In Zoom Out Reset image size

4.2.3. Effects of Increasing the Magnetic Field Strength

Finally, in Figure 9 we investigate the time evolution of models with different initial magnetic-field strengths, B₀ = 0, 1, 3, and 10 μG (initial plasma β₀ = 8.2, 0.92, and 0.082; see Table 2). The panels are the same as in Figure 7, except for an additional panel on the top, showing the Alfvén Mach number $\mathcal {M}_{\rm A}$. Apart from some temporal fluctuations, $\mathcal {M}_{\rm A}$, $\mathcal {M}$, and α_vir are fairly constant over time. Both $\mathcal {M}_{\rm A}$ and $\mathcal {M}$ show some minor systematic decrease, which is caused by dynamo action, amplifying the magnetic field by converting turbulent energy into magnetic energy (Brandenburg & Subramanian 2005). Most of the dynamo action, however, took place already during the first two turbulent crossing times, t < 0 t_ff, during which the turbulence becomes fully established (compare Columns 9 and 14 of Table 2). The dynamo is nearly saturated at t = 0 with only very slow linear amplification happening afterward. In addition, field lines are compressed during local collapse, amplifying the field further in dense cores and clusters (see Figure 5).

Zoom In Zoom Out Reset image size

Most importantly, the last panel of Figure 9 shows that the SFR_ff decreases monotonically with increasing magnetic field because of the stabilizing effect of the magnetic pressure. The strongest magnetic-field case studied here (B₀ = 10 μG, $\mathcal {M}_{\rm A}\approx 1.3$) has an SFR_ff ≈ 0.24, which is almost a factor of two smaller than in the respective purely hydrodynamical run (B₀ = 0, SFR_ff ≈ 0.46). A similar reduction of the SFR with strong magnetic fields compared to purely hydrodynamical or weakly magnetized models is reported in Padoan & Nordlund (2011) and Padoan et al. (2012), who find a maximum reduction by a factor of ∼3. This is a significant but relatively small effect compared to the influence of different forcing on the SFR_ff (see above). Magnetic fields reduce SFR_ff but are unlikely the major player in controlling the SFR, provided that molecular cloud turbulence is super-Alfvénic or at most trans-Alfvénic. This seems to be the case in most clouds. However, as pointed out earlier, on larger scales than molecular clouds, i.e., in the warmer, mainly atomic part of the ISM, turbulence may be trans-Alfvénic or even sub-Alfvénic (Heiles & Troland 2005; Li & Henning 2011; Heyer & Brunt 2012), rendering magnetic fields potentially more important in the process of molecular cloud formation. Still, even inside molecular clouds, magnetic fields seem to reduce fragmentation significantly (see Figure 4), thus potentially having a strong impact on the mass distribution of cores and stars (see also Price & Bate 2007; Hennebelle & Teyssier 2008; Bürzle et al. 2011; Peters et al. 2011; Hennebelle et al. 2011).

We measure Σ_gas and Σ_SFR with a method as close as possible to what observers do to infer Σ_SFR-to-Σ_gas relations (see, e.g., Bigiel et al. 2008; Heiderman et al. 2010), including the effects of telescope beam smoothing. For each simulation, we construct two-dimensional projections of the gas column density Σ_gas and the sink particle column density Σ_SF along each coordinate axis: x, y, z. All maps were smoothed to a resolution N_res/8 with the numerical resolution N_res given in Table 2, such that the size of each pixel in the smoothed maps slightly exceeds the sink particle diameter (which is 5 grid cells; see Section 3.3). We also test smoothing to N_res/32 below, which yields similar results. We then search for pixels with a sink particle column density greater than zero, Σ_SF > 0, and extract the corresponding pixel in the gas column density map, which gives Σ_gas in units of M_☉ pc⁻² for that pixel. The SFR column density is computed by taking the sink particle column density Σ_SF of the same pixel and dividing it by a characteristic timescale for star formation, t_SF, which yields Σ_SFR = Σ_SF/t_SF in units of M_☉ yr⁻¹ kpc⁻². The simplest choice for t_SF is a fixed star formation time, t_SF = 2 Myr, based on an estimate of the elapsed time between star formation and the end of the Class II phase (e.g., Evans et al. 2009; Covey et al. 2010). This is also the t_SF adopted by Lada et al. (2010) and Heiderman et al. (2010) to convert young stellar object (YSO) counts into an SFR column density, so we use it here as the standard approach. However, we also experimented with two other choices for t_SF and present a comparison of those choices below, all yielding similar results.

The result of the procedure explained above is plotted in panel (a) of Figure 11. It shows a scatter plot of Σ_SFR versus Σ_gas measured in all the maps produced from our simulations listed in Table 2 (except for the two low-resolution, 128³ simulations) for a star formation efficiency SFE = 1% (blue) and SFE = 10% (red). Thus, each pixel shown in panel (a) of Figure 11 is one pair of (Σ_gas, Σ_SFR) extracted for each simulation and each projection direction. By combining all data of maps from the three principal projections in x, y, and z, we increase the statistical sample for each model by about a factor of three on average. A total number of 3.5 × 10³ and 1.2 × 10⁴ simulation pixels for SFE = 1% and SFE = 10%, respectively, contribute to the scatter plots in Figure 11. We also add contours of the (Σ_gas, Σ_SFR) distribution, with two contour levels for SFE = 1% (blue contours) and SFE = 10% (red contours). The thick contours enclose 50% of all simulation pixels and the thin contours enclose 99%. The contours help to easily identify the underlying probability distribution of the scattered data points.

Zoom In Zoom Out Reset image size

The simulation data have a broad probability distribution with a clear positive correlation between Σ_SFR and Σ_gas. The data for SFE = 10% are shifted to higher Σ_SFR and lower Σ_gas compared to the SFE = 1% distribution because more gas is accreted by sink particles and thus removed from the gas phase at higher SFE. If we were to fit power laws to the distributions, the slopes would be in the range of 1–2, i.e., Σ_SFR∝Σ^1–2_gas with somewhat flatter slopes at higher SFE.

To compare the simulation data with observations, we add Galactic cloud data from Heiderman et al. (2010) in panel (b) of Figure 11, superimposed on the simulation contours. The observational data are from Galactic observations of clouds and YSOs identified in the Spitzer cores-to-disks (c2d) and GB surveys (Evans et al. 2009) of massive dense clumps (Wu et al. 2010) and of the Taurus molecular cloud (Pineda et al. 2010; Rebull et al. 2010). The simulation data indicated by the same contours of panel (a) fall in the range of the observational data; however, the simulation data show slightly higher Σ_SFR than the observational data, on average. This is not surprising, given that our simulations did not include any local feedback from YSOs. It is known, however, that young stars eject a significant amount of accreted material, thereby reducing the overall accretion rate due to feedback from jets, winds, and outflows (Wardle & Koenigl 1993; Konigl & Pudritz 2000; Beuther et al. 2002; Pudritz et al. 2007; Peters et al. 2011; Seifried et al. 2011a). Hence, only a fraction < 1 of the infalling gas actually ends up on the protostar.

The local core-formation efficiency is parameterized by the factor in Equation (7), from which all the SFR_ff models in Section 2 were derived. Since there is no feedback in our simulations, = 1 by definition. However, we can devise a correction to account for < 1. For this, we simply have to multiply the original Σ_SFR for = 1 by a given < 1. To conserve mass, we also have to account for the fact that a fraction (1 − ) was not accreted and remained in the gas phase due to local feedback. This means we have to increase Σ_gas according to the reduction of Σ_SFR, such that Σ_tot = Σ_gas + Σ_SF with Σ_SF = Σ_SFRt_SF is conserved. Given our simulation data Σ_gas and Σ_SFR with = 1, we can compute values Σ'_SFR and Σ'_gas for < 1 according to the following equations:

$$\begin{eqnarray} \Sigma _{\rm SFR}^\prime (\epsilon) & = & \epsilon \,\Sigma _{\rm SFR}\,, \nonumber \\ \Sigma _{\rm gas}^\prime (\epsilon) & = & \Sigma _{\rm gas}+ \left(1-\epsilon \right) \Sigma _{\rm SF}. \end{eqnarray} \tag{ 47 }$$

Using these expressions, we can correct our simulation data to follow more realistic values of the local efficiency (see also the discussion of in Section 2.3).

The Heiderman et al. (2010) sample of SFR column densities for Galactic clouds shown in panel (b) of Figure 11 is rather broad and presumably covers different evolutionary stages of the clouds, such that a single SFE for the whole sample is quite unlikely. However, since we are currently lacking additional information about the SFE in the observational sample, we can reasonably assume SFEs in the range 1%–10% in the observational data (Evans et al. 2009; Federrath & Klessen 2012, Paper II). In order to find the best-fit local efficiency parameter , we fit our simulated distribution p_sim(Σ'_gas, Σ'_SFR) to the observed distribution p_obs(Σ_gas, Σ_SFR), by applying Equations (47). To do this, we compute the sum of the squared differences Δ² between the two distributions, which have both been sampled to the same (Σ_gas, Σ_SFR) grid with indexes i,

$$\begin{equation} \Delta ^2=\sum _i \big[p_{\rm sim}\big(\Sigma _{\rm gas}^\prime {_{,i}},\,\Sigma _{\rm SFR}^\prime {_{,i}}\big) - p_{\rm obs}\big(\Sigma _{\rm gas}{_{,i}},\,\Sigma _{\rm SFR}{_{,i}}\big)\big]^2,\, \end{equation} \tag{ 48 }$$

for SFE = 1%, 3%, and 10% and for 21 local efficiencies, = [0, 1] in steps of  d = 0.05. For each given SFE, we search for the minimum of Δ² as a function of . This procedure yields best-fit values of the local efficiency parameter = 0.7, 0.5, and 0.3 for SFE = 1%, 3%, and 10%, respectively, in our comparison of simulation data with the Heiderman et al. (2010) Galactic cloud sample.

The simulation data modified to a local efficiency of = 0.5 are shown in panel (c) of Figure 11 together with the original Heiderman et al. (2010) data. Assuming that the observational data have an SFE between 1% and 10%, the local efficiency parameters would be in between = 0.3 and 0.7. This is in good agreement with theoretical models for (Matzner & McKee 2000), with numerical simulations including outflow feedback (Wang et al. 2010; Seifried et al. 2012), and with observational estimates (Beuther et al. 2002, and the discussion on in Section 2.3).

We note that the simulation data in Figure 11 are furthermore consistent with the Galactic cloud samples in Lada et al. (2010) and Gutermuth et al. (2011), showing that Σ_SFR can vary by more than an order of magnitude at any given Σ_gas.

Considering the uncertainties in the SFE from the observations and the uncertainties in the simulations, the overall agreement is encouraging. The HCN(1–0) observational data points of molecular clumps are at the lower end of the distribution, but are still consistent with the simulation data. Possibly, the molecular clumps have a systematically smaller SFE because they are larger structures compared to the YSOs, such that the molecular clumps fall slightly below the general trend. However, this can only be tested when estimates of the cloud SFEs become available (see Paper II, Federrath & Klessen 2012). The Taurus data point as well as a few of the YSO data in the range log₁₀Σ_gas ≈ 1.4–2.8 also lie at the low-Σ_SFR end of the distributions obtained in the simulations. This might be caused by an enhanced magnetic-field influence for these objects. For instance, Taurus seems to be trans-Alfvénic rather than super-Alfvénic (Heyer & Brunt 2012), leading to a reduced Σ_SFR as discussed in Section 4.1. Only one of our MHD simulations approaches this strongly magnetized regime (GT256mM10B10 with $\mathcal {M}_{\rm A}\approx 1.3$; see Table 2), where anisotropies induced by the magnetic field start to become important.

We test the effects of telescope beam smoothing in panel (d) of Figure 11. Panel (d) is identical to panel (c), except that the simulation data were smoothed to grids with resolution N_res/32, i.e., four-times-coarser resolution compared to the contours shown in panel (c). The increased beam smoothing results in distributions with somewhat smaller Σ_gas and Σ_SFR, best seen by comparing the positions of the thickest contours between panels (c) and (d). However, the overall agreement of the simulation data with the Galactic cloud sample is still good, even when the resolution is decreased by a factor of four.

In Figure 12, we study the influence of different choices for the star formation timescale t_SF. The two panels are identical to panel (c) in Figure 11, except for the method by which Σ_SFR = Σ_SF/t_SF was computed in the simulations. The left panel adopts t_SF = t_ff(ρ₀), i.e., the sink particle column density Σ_SF is divided by the freefall time at the mean density ρ₀ of the simulation in which the Σ_SF pixel was found. In the right panel, we use $t_{\rm SF}=t_{\rm ff}(\Sigma _{\rm gas})=\sqrt{3\pi L /(32G\Sigma _{\rm gas})}$, i.e., instead of taking the global mean freefall time, we take the local freefall time of the gas in each pixel. The contours differ slightly between those two last choices and between our standard choice of fixed t_SF = 2 Myr in Figure 11, but the overall agreement between simulation data and Galactic observations is similar in all three cases.

Zoom In Zoom Out Reset image size

Figures 11 and 12 indicate some power-law correlation of the form Σ_SFR∝Σ^N_gas (albeit with significant scatter), similar in exponents N ≈ 1–2 to the Kennicutt–Schmidt relation (Schmidt 1959; Kennicutt 1998) and follow-up measurements for molecular gas (e.g., Wong & Blitz 2002; Gao & Solomon 2004; Bigiel et al. 2008; Kennicutt & Evans 2012). However, the measured values of Σ_SFR in our numerical sample are larger than the extragalactic values of Σ_SFR and larger than theoretical estimates for that regime (e.g., Krumholz et al. 2009) by about one to two orders of magnitude. The Galactic measurements of Σ_SFR by Heiderman et al. (2010) in Figure 11 and by Lada et al. (2010), however, also show values of Σ_SFR that are one to two orders of magnitude above the extragalactic measurements with a scatter of about one to two orders of magnitude. Heiderman et al. (2010) explain this difference between Galactic and extragalactic measurements of Σ_SFR with the different telescope resolutions available for both regimes and thus the different areas over which the measurements of Σ_gas and Σ_SFR are averaged. Both disk-averaged and spatially resolved extragalactic measurements only provide highly smoothed images, mixing both star-forming and non-star-forming gas. Taking these factors into account and correcting for them, Heiderman et al. (2010) conclude that the extragalactic (Σ_gas, Σ_SFR) relations are in agreement with the Galactic measurements. Indeed, decreased telescope resolution (or equivalently observing a region at greater distance) reduces Σ_SFR, but also Σ_gas, as demonstrated here by comparing panels (c) and (d) of Figure 11. Krumholz et al. (2012, p. 71) argue that both Galactic and extragalactic measurements are consistent with a local star formation law, correlating Σ_SFR with Σ_gas/t_SF, where t_SF "is the freefall time evaluated at the density averaged over length scales comparable to the outer scale of turbulence, regardless of the mean density of the region being observed." This seems to be a rather especial definition. Our experiments with three different definitions of t_SF in Figures 11 and 12 do not exclude or prefer any particular choice for t_SF in the Galactic cloud sample studied here. After acceptance of this work, we also learned about a recently submitted paper on a theoretical model for the Σ_SFR–Σ_gas relation by Renaud et al. (2012), which is consistent with our findings for Galactic clouds favoring non-universal behavior of the star formation relation.

The simulations and the observational data shown in Figure 11 are generally in very good agreement. The variations of the observed Σ_SFR in different clouds by up to two orders of magnitude for a given value of Σ_gas (Mooney & Solomon 1988; Lada et al. 2010; Heiderman et al. 2010) and the different scaling relations of Σ_SFR versus Σ_gas (Suzuki et al. 2010) might thus be a result of different physical conditions in Galactic as well as extragalactic molecular clouds. As shown above, star formation is primarily controlled by the forcing and the sonic Mach number of the turbulence, with the magnetic field having a secondary effect. Molecular clouds cover a range of values for these physical parameters and different combinations of those, providing an explanation for the observed scatter in SFRs.

7.1.1. Non-log-normal Effects in the Density PDF

7.1.2. Anisotropies in Sub-Alfvénic Turbulence

7.1.3. Virial Parameter

7.2.1. Approximation of SFR_ff as Constant Over Time

7.2.2. Limited Numerical Resolution

7.2.3. Periodic Boundary Conditions

Assuming a uniform SFE = 1%–10% in the observed Galactic cloud sample by Heiderman et al. (2010), we estimated the local core-formation efficiency parameter = 0.3–0.7 with the best-fit value ≈ 0.5, by fitting our numerical simulations to the observed distribution in Figure 11. There are three major uncertainties in this comparison of the simulations and observations.

First, the SFE in the observed sample is not known. We reasonably assumed SFE = 1%–10%, but some of the individual clouds may not fall in this range. Moreover, there could be a systematic correlation of SFE with gas column density Σ_gas, which is not accounted for. For instance, the HCN(1–0) molecular clump data shown in Figure 11 has potentially smaller SFE on average than the YSO data, because smaller scales tend to exhibit higher SFE (McKee & Ostriker 2007). For instance, it seems plausible that SFE approaches the local core-efficiency , once scales as small as a single core are considered. In contrast, giant molecular cloud complexes as a whole typically only have SFEs of a few percent at most (see Paper II, Federrath & Klessen 2012).

The second uncertainty is the effect of the telescope resolution. Lower resolution (or observation of a very distant region, e.g., a whole galaxy) inevitably means that the observed star-forming regions are smoothed over larger areas compared to a high-resolution observation of the same region. The effect of reducing the beam resolution by a factor of four in our synthetic observations of the simulated clouds is demonstrated by comparing panels (c) and (d) in Figure 11, resulting in a relatively weak but noticeable reduction of Σ_SFR and Σ_gas.

The third major uncertainty is the star formation timescale t_SF used to convert a given star formation column density Σ_SF into a rate Σ_SFR = Σ_SF/t_SF. In Figure 11, we adopted a fixed t_SF = 2 Myr as often used by observers (e.g., Heiderman et al. 2010; Lada et al. 2010). However, we studied two additional choices of t_SF in Figure 12: one where t_SF = t_ff(ρ₀) (division by the global freefall time) and the other where t_SF = t_ff(Σ_gas) (division by the local freefall time). Comparing these three choices for t_SF, we find that the resulting Σ_SFR-to-Σ_gas correlations change slightly, but the overall effect is rather weak. Given the broad distributions in both the simulation data and in the Heiderman et al. (2010) Galactic cloud sample, it is hard to decide which method provides better agreement. They all seem to agree reasonably well within the observational range of Galactic clouds.

Finally, we note a fundamental difficulty of estimating actual SFRs or SFR_ff in observations. Cloud observations are inevitably limited to a nearly instantaneous snapshot of the state of a cloud with respect to the relevant timescales for star formation, which exceed the lifetime of a human being by orders of magnitude. However, measuring a real SFR requires knowledge about the time evolution of the cloud, which is thus not available. Strictly speaking, a direct measurement of the time derivative of star formation, i.e., the SFR is thus impossible in observations. This is why we can only make meaningful comparisons of star formation in simulations and observations based on the methods explained and applied in Section 6 (Figures 11 and 12), but not the actual SFRs computed from the time evolution of star formation.

• 1.  
  In Section 2, we derived six analytic models for the SFR per freefall time, SFR_ff: the original Krumholz & McKee (2005, KM), Padoan & Nordlund (2011, PN), and Hennebelle & Chabrier (2011, HC) models and the multi-freefall KM, PN, and HC models, which are all based on an integral over the density PDF, Equation (1), leading to different analytic solutions for SFR_ff, summarized in Table 1. They all yield a dimensionless SFR per freefall time, SFR_ff, based on Equation (7), which can be transformed to a real SFR with units of M_☉ yr⁻¹ by applying Equation (6).
• 2.  
  We extended the (multi-freefall) KM and (multi-freefall) HC theories to include magnetic fields by introducing a magnetic-pressure correction given by Equation (17), which allows us to replace the sound speed by an effective magnetosonic speed given by Equation (18) or (19) for super-Alfvénic, isothermal turbulence.
• 3.  
  We analyzed the basic dependences of all six theories on the four parameters listed above. SFR_ff decreases with increasing virial parameter α_vir, while it increases with increasing sonic Mach number $\mathcal {M}$ in the best multi-freefall theories (see Figure 1). Varying the forcing parameter b from purely solenoidal forcing (b = 1/3) to purely compressive forcing (b = 1) leads to a higher SFR_ff by more than an order of magnitude (Figure 2). Stronger magnetic fields parameterized by decreasing plasma β (or equivalently decreasing Alfvén Mach number $\mathcal {M}_{\rm A}$) lead to decreasing SFR_ff (Figure 3).

• 1.  
  In Sections 3 and 4, we performed a set of numerical experiments of star formation, covering molecular cloud sizes and masses in the range L = 0.3 to 200 pc and M_c = 300 to 4 × 10⁶ M_☉ (see Table 2) with solenoidal, mixed, and compressive forcings of the turbulence (see Section 3.2 for details of the forcing) to test the analytic models. We also ran super-Alfvénic simulations with varying magnetic-field strength to test the influence of magnetic fields on the SFR. All simulations include sink particles to model core and star formation, allowing us to measure SFR_ff, depending on α_vir, $\mathcal {M}$, b, and β.
• 2.  
  We computed the virial parameter α_vir ≡ 2E_kin/|E_grav| based on the uniform-density, spherical approximation given by Equation (15), and based on the actual, three-dimensional, inhomogeneous gas distribution in the simulations. Depending on the forcing and Mach number of the turbulence, we find that these two definitions can differ by more than an order of magnitude (compare Columns 10 and 11 in Table 2), which means that theoretical and observational estimates of α_vir based on a uniform-density, spherical approximation must be considered with caution.
• 3.  
  The SFR converges with increasing numerical resolution (Figures 7 and 10). The statistical uncertainty in SFR_ff is about a factor of 1.5, indicated by comparing three different random realizations of the same parameter set (Figure 7), similar to the uncertainty introduced by limited numerical resolution.
• 4.  
  We found that for our models with $\mathcal {M}\sim 10$, compressive forcing yields SFRs at least an order of magnitude higher than solenoidal forcing, emphasizing the role of different turbulent energy injection mechanisms for the SFR (Figure 7). The cloud morphology also depends strongly on the type of forcing and sonic Mach number (see Figure 6). The SFR increases by a factor of about four for compressive forcing between $\mathcal {M}=3$ and $\mathcal {M}=50$ (Figure 8).
• 5.  
  Including magnetic fields in simulations with $\mathcal {M}\sim 10$ and mixed turbulent forcing, we found that the magnetic field is amplified in regions of core and cluster formation (Figure 5), reducing the SFR_ff by a factor of two between purely hydrodynamic turbulence ($\mathcal {M}_{\rm A}\rightarrow \infty$) and trans-Alfvénic turbulence with $\mathcal {M}_{\rm A}\sim 1.3$ (see Figure 9). This is a relatively small change in SFR_ff for such a fairly strong magnetic field, compared to the dependence of SFR_ff on α_vir, $\mathcal {M}$, and b. However, magnetic fields do affect the morphology of the clouds even on large scales, and they reduce fragmentation (see Figure 4), thus potentially having an important impact on the core and stellar IMFs.
• 6.  
  A detailed comparison of SFR_ff (simulation) with SFR_ff (theory) in Figure 10 showed that the multi-freefall analytic theories are generally better than the non-multi-freefall theories. The multi-ff KM and multi-ff PN models give the best fits to our simulation data (see Tables 3 and 4) with reasonable best-fit model parameters, 1/ϕ_t ≈ 0.5 for both multi-ff KM and PN models, as well as ϕ_x ≈ 0.17 for the multi-ff KM model, and θ ≈ 1 for the multi-ff PN model, suggesting a close connection between the magnetothermal Jeans scale and the magnetosonic scale, as well as turbulence driven on the outer, largest scales of molecular clouds.
• 7.  
  All numerical simulations agree with the multi-ff KM and PN theories within a factor of three, and come closer to the analytic prediction with increasing numerical resolution. This is an encouraging agreement, given that the modeled SFRs vary over two orders of magnitude in our numerical simulations (see Figure 10).

• 1.  
  We compared our numerical simulations with observations of the SFR column density Σ_SFR as a function of the gas column density Σ_gas, measured in Galactic clouds in Section 6 (Figure 11). We showed that the simulations slightly overestimate the SFR compared to the observed clouds because we did not include any local radiative and mechanical feedback from young stellar objects, and hence, the local efficiency parameter = 1 in our simulations, by definition. However, assuming a constant, global star formation efficiency in the observed clouds of SFE ≈ 1%–10% (see Paper II, Federrath & Klessen 2012), we can adjust our numerical simulation data with Equations (47) to account for < 1. Doing so, we found the best-fit local efficiency ≈ 0.5 ( = 0.7, 0.5, and 0.3 for SFE = 1%, 3%, and 10%, respectively) for the observed Galactic clouds, which is in good agreement with theoretical expectations, independent numerical simulations, and observations of individual protostellar cores.
• 2.  
  We studied the effects of telescope beam smoothing in panels (c) and (d) of Figure 11, and the effect of varying the definition of the star formation timescale t_SF to determine Σ_SFR in Figure 12. We found that both the telescope beam resolution and the definition of t_SF introduce minor uncertainties in our comparison between simulations and observations.
• 3.  
  The correlation between Σ_gas and Σ_SFR in Figure 11 is consistent with power laws of the form Σ_SFR∝Σ^N_gas with exponents N = 1–2 (albeit with significant scatter), which is similar to extragalactic measurements of Σ_gas–Σ_SFR correlations.

The overall agreement between theory, simulations, and observations in Figures 10 and 11 is encouraging, considering the simplifications inherent in the theoretical models, the limitations of the numerical simulations, and the uncertainties in the SFEs of the observed sample of clouds (see Section 7). We conclude that supersonic, magnetized turbulence is a key process, likely controlling the SFR of molecular clouds in the Milky Way and potentially in other galaxies.