THE DENSITY VARIANCE–MACH NUMBER RELATION IN SUPERSONIC, ISOTHERMAL TURBULENCE

The last few years have seen an increasing number of analytic models of the star formation process that use the log-normal density probability distribution function (PDF) produced by supersonic turbulent flows to predict statistical quantities such as the initial and/or core mass function (e.g., Padoan & Nordlund 2002; Hennebelle & Chabrier 2008, 2009) and the star formation rate (Krumholz & McKee 2005; Padoan & Nordlund 2009). A key assumption in these models is a relationship identified in early numerical studies between the PDF width—the density variance or standard deviation—and the root-mean-square (rms) Mach number $\mathcal {M}$ in supersonic, isothermal turbulence. The relationship is generally assumed to be linear in the standard deviation of linear density, i.e.,

$$\begin{equation} \sigma _{\rho /\bar{\rho }} =b\mathcal {M}, \end{equation} \tag{ 1 }$$

where b is a constant of order unity and density is scaled in terms of the mean, $\bar{\rho }$. For a log-normal distribution, this is equivalent to

$$\begin{equation} \sigma _{s}^{2}=\ln \left(1+b^{2}\mathcal {M}^{2}\right), \end{equation} \tag{ 2 }$$

where $s\equiv \ln (\rho /\bar{\rho })$, such that σ_s is the standard deviation in the logarithm of density.

Apart from the early empirical findings of Vazquez-Semadeni (1994), Padoan et al. (1997b), and Passot & Vázquez-Semadeni (1998), there is no clear reason why the relationship should be of this form. Mathematically, the appearance of a log-normal distribution can be understood as a consequence of the multiplicative central limit theorem assuming that individual density perturbations are independent and random (Vazquez-Semadeni 1994; Passot & Vázquez-Semadeni 1998; Nordlund & Padoan 1999). In physical terms this has been interpreted as meaning that density fluctuations at a given location are constructed by successive passages of shocks with a jump amplitude independent of the local density (e.g., Ballesteros-Paredes et al. 2007; Kritsuk et al. 2007; Federrath et al. 2010). However, it has not so far proved possible to analytically predict the relationship based on these ideas (though see Padoan & Nordlund 2009). Thus, a common approach in numerical studies of turbulence—usually at a fixed Mach number—has been to measure the parameter b, assuming Equation (2), that gives best-fitting log-normal to the time-averaged PDF. However, reported estimates for b are widely discrepant. For example, Padoan et al. (1997b) found b ≈ 0.5 while more recently Kritsuk et al. (2007) (at Mach 6) find a much lower value of b ≈ 0.26 and Beetz et al. (2008) find b ≈ 0.37, while Passot & Vázquez-Semadeni (1998) found b ≈ 1 (though with some confusion over σ_s versus $\sigma _{\rho /\bar{\rho }}$).

Federrath et al. (2008) and Federrath et al. (2010) reconcile these results in part by the finding that the width of the PDF depends not only on the rms Mach number but also on the relative degree of compressible and solenoidal modes in the turbulence forcing, with b = 1/3 appropriate for purely solenoidal and b = 1 for purely compressive forcing. This is in keeping with earlier discussions by Passot & Vázquez-Semadeni (1998) and Nordlund & Padoan (1999, p. 222), the latter authors noting that "for compressional forcing at low Mach numbers (leading to an ensemble of sound waves), the standard deviation is expected to be equal to the rms Mach number itself."

Observationally, the log-normality of the three-dimensional (3D) PDF is reflected in the two-dimensional (2D) column density PDF, for example, as measured from dust extinction maps (e.g., Lombardi et al. 2006, 2008, 2010; Kainulainen et al. 2009)—at least in earlier stages of molecular cloud evolution, suggesting that this phase could be dominated by roughly isothermal turbulence in which self-gravity is relatively unimportant. Only for seemingly more evolved clouds (including Taurus) do Kainulainen et al. (2009) see significant tails at higher (column) densities (similarly found by Lombardi et al. 2010). However, measurements of the projected 2D variance (or PDF) cannot be directly used to constrain the relationship with Mach number. Recently, Brunt et al. (2010a, hereafter BFP) have shown how projection effects can be overcome to infer the 3D density variance from column density observations, in turn leading to a method for extracting the unprojected (3D) density PDF from the observational data (Brunt et al. 2010b). This enables the relationship between the standard deviation in linear density and Mach number to be tested observationally, with initial application to Taurus finding b = 0.48^+0.15_−0.11 (Brunt 2010). A similar method was employed by Padoan et al. (1997a) to infer the 3D density variance from extinction measurements in IC5146, similarly finding b ≈ 0.5.

The problem with all of the above is that calculations—or observations—performed at a single (rms) Mach number can only ever assume the relationship given by Equation (2) or Equation (1) and cannot be used to constrain it unless a range of Mach numbers are studied. Indeed, Lemaster & Stone (2008,, hereafter LS08)—performing a series of calculations with Mach numbers in the range $1.2\le \mathcal {M}\le 6.8$—find a relationship

$$\begin{equation} \sigma ^{2}_{s}=0.72\ln \left(1+0.5\mathcal {M}^{2}\right)-0.20, \end{equation} \tag{ 3 }$$

based on a fit to measurements of the mean in the logarithm of density, $\bar{s}$, as a function of $\mathcal {M}$, which we have here converted to a σ_s–$\mathcal {M}$ relation using $\bar{s}=-\sigma _{s}^{2}/2$. However, this three-parameter fit is clearly not unique, and it remains to be determined whether this or a similar relationship continues to hold at higher Mach numbers.

The present study is motivated by a need to compare the theoretical predictions with the observational constraints. In particular, LS08 only perform calculations up to $\mathcal {M}\approx 6.8$—corresponding to one-dimensional (1D) line full width at half-maximum of ∼1.6 km s⁻¹ (at 10 K), which is rather low in terms of what is found in the real interstellar medium (ISM). In particular, Taurus has $\mathcal {M}\sim 17$, so a study going up to (at least) Mach 20 or so is needed. Our aim in this Letter is precisely this: to pin down the theoretical relationship—with as few assumptions as possible—up to sufficiently high Mach numbers that a meaningful comparison can be made with observed molecular clouds. While additional physics such as non-isothermality (e.g., Scalo et al. 1998), the multiphase nature of the ISM, and self-gravity (e.g., Klessen 2000; Kritsuk et al. 2010) are all expected to change the theoretical predictions at some level, the isothermal, non-self-gravitating case is an important reference point that remains theoretically uncertain. Furthermore, a clear prediction for this simple case can be used to gauge the relative importance of such additional physics in observed clouds.

2.1. Log-normal Distributions

The log-normal distribution is given by

$$\begin{equation} p(s)ds=\frac{1}{\sqrt{2\pi \sigma _{s}^{2}}}\exp {\left[-\frac{1}{2}\left(\frac{s-\bar{s}}{\sigma _{s}}\right)^{2}\right]}ds, \end{equation} \tag{ 4 }$$

where $s\equiv \ln (\rho /\bar{\rho })$ such that $\bar{s}$ and σ_s denote the mean and standard deviation in the logarithm of (scaled) density, respectively, and $\bar{\rho }$ is the mean in the linear density. The mean and variances in a log-normal distribution are related by

$$\begin{equation} \bar{s}=-\frac{1}{2}\sigma _{s}^{2}, \end{equation} \tag{ 5 }$$

and

$$\begin{equation} \sigma _{\rho /\bar{\rho }}^{2}= \exp \left(\sigma ^{2}_{s}\right)-1. \end{equation} \tag{ 6 }$$

2.2. Numerical Simulations

2.3. Measuring the Density Variance

3.1. σ_s as a Function of $\mathcal {M}$

The direct measurements of the standard deviation in the log density, σ_s, are shown in Figure 1 from the results of calculations performed using 128³, 256³, and 512³ SPH particles, using methods (2) and (4) (see the legend). Dashed lines show the standard relation (Equation (2)) with b = 1/3 and b = 1/2, while the dotted line shows the best-fitting relationship found by LS08 (Equation (3)). Both the b = 1/3 curve, expected for solenoidally driven turbulence (Federrath et al. 2008, 2010) and the LS08 fit show reasonable fits to the data, and indeed cannot be distinguished given the time variability present in the calculations. However, adopting b = 0.5 is clearly not consistent with our solenoidally driven results. The results are also consistent with our earlier findings (PF10), that showed a convergence in both grid and SPH methods toward b ≈ 0.35–0.4 at Mach 10.

Zoom In Zoom Out Reset image size

The measured value of σ_s cannot be compared to observations, since only the 3D variance in the linear density can be observationally inferred (e.g., using the BFP method). Thus, it is necessary to either measure or infer the linear variance from simulations to make this comparison.

3.2. Direct Measurement of $\sigma _{\rho /\bar{\rho }}$

A direct measurement of the linear density variance is difficult even with high-resolution simulations, demonstrated in Figure 2, that shows the directly measured $\sigma _{\rho /\bar{\rho }}^{2}$ as a function of σ²_s. The dashed line shows the expected relationship for a log-normal distribution (Equation (6)). The exponential relationship is resolved only at the lowest Mach numbers, $\mathcal {M}\lesssim 5$, corresponding to σ²_s ≲ 1.4, while at higher Mach numbers $\sigma _{\rho /\bar{\rho }}$ is a strong function of resolution, most easily demonstrated by interpolating the highest resolution SPH calculations (512³ particles) to fixed grids of decreasing resolution (see the legend). This dependence is the reason why we eventually interpolated the SPH data onto an adaptive mesh, refined such that Δx < h for all cells within the smoothing radius 2h of any given particle, in order to obtain a result for $\sigma _{\rho /\bar{\rho }}$ that captures the maximum resolution available in the SPH simulations, though we remain limited by the intrinsic resolution of the simulations. PF10 found that SPH simulations at Mach 10 resolved a maximum density at 128³ particles similar to that captured on a fixed grid at 512³ grid cells. This is consistent with the results here, where it is necessary to refine the grid to an effective 8192³ for the Mach 20 calculations employing 512³ SPH particles. It is also evident that fully resolving the strong fluctuations in the linear density at high Mach number is intractable with current computational resources. Our findings also suggest that the linear density variance is likely to be severely underestimated by limited observational resolution, so the results in Taurus and IC5146 are almost certainly lower limits (see also Brunt 2010).

Zoom In Zoom Out Reset image size

3.3. $\sigma _{\rho /\bar{\rho }}$ as a Function of $\mathcal {M}$: Comparison to Observations

The direct—but resolution limited—measurements for $\sigma _{\rho /\bar{\rho }}$, computed without assumptions from the adaptive mesh refinement (AMR) grid, are shown in Figure 3 (see the legend,i.e., method (1) above). A better approach is to use the fact that the measurements for σ_s are resolution independent (Figure 1), meaning that we can use the assumption of log-normality to infer the fully resolved value for $\sigma _{\rho /\bar{\rho }}$, i.e., using Equation (6) (method 2). The standard deviations computed in this way are also shown and as expected are consistent with a linear $\sigma _{\rho /\bar{\rho }}$–$\mathcal {M}$ relation with b = 1/3 for solenoidal forcing (Federrath et al. 2008, 2010), and also with the LS08 best fit, the latter translated in terms of $\sigma _{\rho /\bar{\rho }}$ using Equation (6).

We are now in a position to compare with the observational results. It should first be noted that the assumption of the σ_s–$\sigma _{\rho /\bar{\rho }}$ relation for the simulations essentially gives us an upper limit on $\sigma _{\rho /\bar{\rho }}$, whereas (see above) the finite resolution of the observations almost certainly gives a lower limit on the variance. The results from Padoan et al. (1997a; in IC5146: b = 0.5 ± 0.05 at Mach 10) and Brunt (2010; in Taurus, $\mathcal {M}_{\rm {rms}}=17.6\pm 1.8$ and $\sigma _{\rho /\bar{\rho }}=8.49^{+1.85}_{-1.71}$) are also plotted, both of which are consistent with b = 0.5 but clearly inconsistent with the calculations of purely solenoidally driven turbulence.

The other extreme is given by the b = 1 line in Figure 3, corresponding to purely compressive forcing (Federrath et al. 2008, 2010). Data points from two 1024³ grid simulations of purely solenoidal and purely compressive forcing by Federrath et al. (2008, 2010) are also shown (cyan triangles). Clearly, purely solenoidal and purely compressive forcing seem inconsistent with the observations, while a mixture and/or the addition of gravity can fit the observations, best fit by a linear relation with b ≈ 0.5. The σ_ρ–$\mathcal {M}$ plane, however, needs to be populated with many more observational measures to draw more definite conclusions about, e.g., regional and evolutionary variations.

Zoom In Zoom Out Reset image size

Since molecular clouds are observed to be magnetized, we have additionally computed a series of magnetohydrodynamics (MHD) calculations, driven identically to their hydrodynamic counterparts but with an initially uniform magnetic field threading the box. These have been computed on a fixed grid using the flash code, as described in Federrath et al. (2010), PF10, and Brunt et al. (2010a, 2010b; note that the driving routines are implemented identically in both the SPH and grid code). The magnetic field strength in the MHD calculations is characterized by the ratio of gas-to-magnetic pressure β = P/P_mag in the initial conditions, where $P_{{\rm mag}}=\frac{1}{2} B^{2}/\mu _{0}$. As done previously, we have also examined the effect of resolution, with a finding similar to the hydrodynamic case, namely, that direct measurements of $\sigma _{\rho /\bar{\rho }}$ are strongly resolution affected (underestimated)—even at modest Mach numbers, similar to the results shown in Figure 2, but that measurements of σ_s are resolution independent (at ≳256³ grid cells).

Figure 4 shows the results, similar to Figure 1 for the hydrodynamic case and overplotted with the hydrodynamic b = 1/3 and b = 1/2 relationships (dashed lines) as well as the best-fitting MHD relationship found by LS08. The most obvious difference with MHD is that the density variances are significantly lower than their hydrodynamic counterparts at high (sonic) Mach number $\mathcal {M}\gtrsim 10$—though conversely marginally higher at lower Mach numbers. On the whole increasing the field strength seems to decrease the mean σ_s slightly. While a complete MHD study is beyond the scope of this Letter, the results clearly illustrate a shallower relationship than the hydrodynamic b = 1/3 curve at high Mach number, with σ_s only weakly dependent on $\mathcal {M}$ in this regime (for β ≲ 1). Thus, if anything, adding magnetic fields decreases the variance in the density field, worsening the discrepancy with observations.

Zoom In Zoom Out Reset image size

We have measured the relationship between the density variance and the mean-square Mach number from a series of simulations of supersonic, isothermal, solenoidally driven turbulence over a wide range of Mach numbers ($1\lesssim \mathcal {M}\lesssim 20$). We find that the standard relationship given by Equation (2) with b = 1/3 provides a good fit to the data over this range, consistent with the heuristic theory of Federrath et al. (2010) for solenoidal driving and with similar measurements by Lemaster & Stone (2008) at lower Mach numbers. While it is difficult to measure the variance in linear density directly from simulations with finite resolution, the inferred relationship (Equation (1), with b = 1/3) appears inconsistent with observational determinations (b ≈ 0.5) in Taurus and IC5146, suggesting that additional physics such as gravity is important in these clouds and/or that some form of compressive driving is relevant. This is consistent with the findings of Kainulainen et al. (2009) and Lombardi et al. (2010) for Taurus, where self-gravity is invoked to explain the deviation from log-normality in the high density tail of the (column) density PDF. Magnetic fields do not help to explain the discrepancy.

We acknowledge computational resources from the Monash Sun Grid (with thanks to Philip Chan) and the Leibniz Rechenzentrum Garching. Plots utilized splash (Price 2007) with the new giza back end by James Wetter. C.F. is supported by Baden-Württemberg Stiftung grant P-LS-SPII/18.