THE COLUMN DENSITY VARIANCE IN TURBULENT INTERSTELLAR MEDIA: A FRACTAL MODEL APPROACH

A variety of observations and simulations have shown that the density structures of the interstellar medium (ISM) are scale-free, hierarchical, and fractal (e.g., Stutzki et al. 1998; Elmegreen & Scalo 2004; Burkhart et al. 2012). The supersonic, compressible turbulence is likely responsible for the complex and hierarchical density structures observed in the ISM. The density fluctuation due to the supersonic turbulence plays a crucial role in models of star formation rates and initial mass functions (Mac Low & Klessen 2004; Krumholz & McKee 2005; Hennebelle & Chabrier 2008; Federrath & Klessen 2012). Hence, the statistical properties of density structures in the turbulent ISM have been extensively studied in hydrodynamic (HD) and magnetohydrodynamic (MHD) simulations.

It is now well known that the probability distribution functions (PDFs) of the three-dimensional (3D) densities and column densities of the turbulent ISM are close to lognormal (Vázquez-Semadeni 1994; Nordlund & Padoan 1999; Klessen 2000; Ostriker et al. 2001; Wada & Norman 2001; Burkhart & Lazarian 2012). The standard deviation of density ($\sigma _{\rho /\rho _{0}}$) increases with sonic Mach number, M_s (Nordlund & Padoan 1999; Ostriker et al. 2001). This trend is expressed by $\sigma _{\rho /\rho _{0}}=bM_{{\rm s}}$, where ρ/ρ₀ is the density (ρ) normalized by the mean density (ρ₀) and b is a constant of proportionality, known as the turbulence forcing parameter. Values of b ranging from ≈0.3 to ≈1.0 have been suggested in numerical simulations (Padoan et al. 1997; Passot & Vázquez-Semadeni 1998; Kritsuk et al. 2007; Beetz et al. 2008). Federrath et al. (2008, 2009, 2010) have found that for the same M_s, compressive (dilatational) forcing leads to much larger density variance compared to solenoidal (rotational) forcing. They showed that b depends on the type of the turbulence forcing varying from ≈1/3 for solenoidal forcing to ≈1 for compressive forcing.

The density power spectrum is also a useful tool to characterize the turbulent ISM. The power spectra of turbulent clouds are power laws in form (Crovisier & Dickey 1983; Despande et al. 2000; Stanimirović & Lazarian 2001; Padoan et al. 2004). Kim & Ryu (2005) have investigated the dependency of the power-law slope of the density power spectrum on M_s and found that the density spectrum becomes gradually shallower as M_s increases in HD turbulent media. Kowal et al. (2007) and Burkhart et al. (2010) confirmed the flattening of the density spectra with M_s in MHD cases. Federrath et al. (2009) showed that compressive forcing leads to significantly steeper density spectra than solenoidal forcing.

We note that the spectral slope of the 3D density power spectrum is the same as that of the projected column density power spectrum (e.g., Stutzki et al. 1998; Padoan et al. 2004) and thus can be directly extracted from observations. On the other hand, the 3D density is not an observationally accessible quantity, but instead the column density, i.e., the integral of the density along a line of sight, is observable. There is a great deal of observational data from surveys, which trace the ISM column density. Some examples are the GALFA (Peek et al. 2011) and GASKAP (Dickey et al. 2012) H i surveys and the Columbia-Cfa (Dame et al. 2001) and ThrUMMS¹ molecular gas surveys. The best way to infer the variance of the 3D density field from observations is obviously to compare the observational data with a large number of turbulence simulations. However, instead of performing expensive simulations, a simpler approach to derive the 3D density variance would be of great interest. Padoan et al. (1997) generated random density distributions with various density variances and spectral slopes and compared the resulting synthesized extinction maps with the observed extinction data. Brunt et al. (2010) developed an interesting technique for calculating the variance of a 3D density field using only information contained in the two-dimensional projection. The technique has been applied to the Taurus molecular cloud (Brunt 2010). More recently, Burkhart & Lazarian (2012) found a column density variance–Mach number relationship, using MHD simulations. The relationship closely follows the form of the 3D density variance–Mach number relationship but includes a scaling parameter A such that

$$\begin{equation} \sigma _{\ln N}^{2}=A\ln \big(1+b^{2}M_{{\rm s}}^{2}\big), \end{equation} \tag{ 1 }$$

where σ²_ln N is the variance of the logarithm of the column density N.²

In the present study, quantitative expressions relating the 3D density variance, spectral slope, and column density variance are obtained by using the fractional Brownian motion (fBm) algorithm (Saupe 1988; Stutzki et al. 1998; Elmegreen 2002). A lognormal density field is obtained by generating a Gaussian random field with a power-law spectral slope (γ_g) and exponentiating the field (Elmegreen 2002). This Letter is organized as follows. In Section 2, we describe how to infer the variance of the 3D density field from the given spectral slope and column density variance. Section 3 investigates the σ²_ln N–M_s relationship. Concluding remarks are given in Section 4.

The fBm structures are generated by first assigning 3D Fourier coefficients following the prescription of Elmegreen (2002). The Fourier amplitudes are generated to be distributed as a normal Gaussian with the variance $\langle \left|a(\mathbf {k})\right|^{2}\rangle \ =\ \mathbf {|k|}^{-\gamma _{{\rm g}}}$. The inverse Fourier transform gives a Gaussian random field ρ_g($\mathbf {x}$). We then multiply the density field with the desired standard deviation of the logarithmic density σ_ln ρ, and exponentiate the field to obtain ρ = exp (σ_ln ρρ_g).

The standard deviation σ_ln N depends on both γ_g and σ_ln ρ. We should also note that the spectral slopes γ of resulting lognormal density fields are different from the input slopes γ_g. Therefore, we need to find the relationship between the input variables (γ_g, σ_ln ρ) and the outputs (γ, σ_ln N). In addition, the density fields generated with the same parameters (γ_g, σ_ln ρ) show large fluctuations due to random phases and amplitudes in Fourier space. We thus generated a large volume of lognormal density fields by varying γ_g and σ_ln ρ and obtained the average relationships between the parameters. We varied γ_g from 2.0 to 5.25 in steps of 0.25 and σ_ln ρ from 0.01 to 2.76 in steps of 0.25. One hundred random realizations with a box size of 128³ were generated for each combination of γ_g and σ_ln ρ, resulting in a total of 16,800 realizations. We also generated 30 realizations with a resolution of 256³ for every (γ_g, σ_ln ρ) pair and confirmed that the 256³ data set is consistent with the results obtained with the 128³ resolution. As will be shown later, the present results are also consistent with the limited samples obtained with a box size of 512³. The standard deviations and spectral slopes were calculated for every realization and averaged to obtain the ensemble average values of (γ, σ_ln N) for each combination of the input (γ_g, σ_ln ρ) values. The power spectra of the resulting lognormal density fields show a slight flattening at large wavenumbers. Examples of the spectra are shown in Figure 1. Spectral slopes were obtained by least-squares fits over a wavenumber range of 3 ⩽ k ⩽ 35. Here, the dimensionless wavenumber is defined by k = L/λ with the wavelength λ and box size L.

Zoom In Zoom Out Reset image size

For ease of use, we parameterized the relationships between the parameters γ, γ_g, σ_ln ρ, and σ_ln N using polynomials. For every σ_ln ρ, σ_ln N was parameterized with a cubic function of γ, i.e.,

$$\begin{equation} \sigma _{\ln N}=\sum _{i=0}^{3}c_{i}\gamma ^{i}=c_{0}+c_{1}\gamma +c_{2}\gamma ^{2}+c_{3}\gamma ^{3}. \end{equation} \tag{ 2 }$$

The coefficients c_i (i = 0, ..., 3) were then fitted with a quartic polynomial of σ_ln ρ, i.e.,

$$\begin{equation} c_{0}=\sum _{j=0}^{4}p_{j}\sigma _{\ln \rho }^{j}=p_{0}+p_{1}\sigma _{\ln \rho }+p_{2}\sigma _{\ln \rho }^{2}+p_{3}\sigma _{\ln \rho }^{3}+p_{4}\sigma _{\ln \rho }^{4}. \end{equation} \tag{ 3 }$$

In this way, σ_ln N could be estimated as a function of γ and σ_ln ρ. We also parameterized γ with a cubic function of γ_g and γ_g with another cubic function of γ for every σ_ln ρ, i.e., γ = Σ_ia_iγⁱ_g and γ_g = Σ_ib_iγⁱ (i = 0, ..., 3). The coefficients a_i and b_i (i = 0, ..., 3) were then fitted with quartic polynomial functions of σ_ln ρ, as done for the coefficients c_i.

The final coefficients are shown in Table 1. Figure 2 shows the resulting γs as functions of γ_g for every σ_ln ρ. Figure 3 presents the resulting (γ, σ_ln N) pairs for every (γ_g, σ_ln ρ) pair. In the figures, the solid and dotted lines are the best-fit polynomials to reproduce the (γ, σ_ln N) pairs corresponding to constant σ_ln ρ and γ_g, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The power spectra of lognormal density (ρ) fields are usually shallower than those of the original Gaussian density (ln ρ) fields, i.e., γ < γ_g. The property that the power spectra of ln ρ are stiffer than those of ρ in turbulence simulations is also noticeable in Table 2 of Kowal et al. (2007). The difference between γ and γ_g becomes more significant for a larger σ_ln ρ and for a shallower (smaller) γ_g. Column density variance is always smaller than the variance of 3D density, as expected. The ratio of σ_ln N to σ_ln ρ becomes smaller for a shallower γ. This is caused by less spatial correlation in the density fields with a shallower spectral slope. For instance, the density fields with γ = 0 have no spatial correlation.

The equations relating γ and γ_g will be useful in generating the realistic ISM density structures. The equation relating σ_ln N to γ and σ_ln ρ or Figure 3 provides a tool to derive σ_ln ρ from the observable σ_ln N and γ. Once we obtain γ and σ_ln N from observations, we can vary σ_ln ρ to find the best-fit value that gives the observed σ_ln N.

We now investigate the relationship between σ_ln N and M_s by using the expressions derived in Section 2. In the above, we have assumed that γ and σ_ln ρ are independent variables. However, γ and σ_ln ρ are not independent in the real turbulent media, but rather depend on M_s. If we ignore magnetic fields, the statistical properties are uniquely determined by M_s and b. We therefore need to specify the dependencies of γ and σ_ln ρ on M_s.

Seon (2009) combined the results of Kim & Ryu (2005), Kritsuk et al. (2006), and Padoan et al. (2004) and obtained a simple relation between M_s and γ, i.e., γ = 3.81 M^−0.16_s, which may be applicable to solenoidal forcing. The power spectra of the compressively driven turbulence are considerably steeper than those of the solenoidal forcing. The spectral slopes at M_s ≈ 2.3 and 5.6 for compressive forcing are provided by Schmidt et al. (2009) and Federrath et al. (2009), respectively, and are ∼0.6 larger than the values for the solenoidal case. We therefore assume that γ = 3.81 M^−0.16_s + 0.6 for the compressive mode. There are two transverse (solenoidal) and one longitudinal (compressive) wave modes in the 3D space. Therefore, the solenoidal and compressive forcing would be mixed with a ratio of 2:1 if no preferential forcing mode were provided (Federrath et al. 2010). The power spectrum of this natural mixing mode might be obtained by adding two power spectra for solenoidal and compressive forcing with a ratio of 2:1. The resulting power spectrum is also represented by a power law with a spectral slope, which is at most ∼0.06 larger than the slope of the solenoidal case. We therefore assume the same spectral slope as the solenoidal forcing case for the natural mixing case.

The density variance in logarithmic scale is obtained from the variance in linear scale using the property of a lognormal distribution, i.e., σ²_ln ρ = ln (1 + σ_ρ/ρ0²) = ln (1 + b²M²_s). The forcing parameter b is 1/3, 0.4, and 1 for solenoidal, natural mixing, and compressive modes, respectively (Federrath et al. 2008, 2009, 2010).

Using γ and σ_ln ρ as functions of M_s, and Equation (2), we calculate σ_ln N for solenoidal, natural mixing, and compressive modes as functions of M_s. The results are shown in Figure 4, wherein the upper three black curves show σ_ln N as functions of M_s for three different forcing types. We fitted the obtained curves with Equation (1). The best-fit A values are found to be 0.2, 0.22, and 0.38 for the solenoidal, natural mixing, and compressive modes, respectively. In Figure 4, the b and best-fit A values are denoted in parentheses along the curves. The best-fit curves are shown as dashed lines.

Zoom In Zoom Out Reset image size

Applying the expression relating γ and γ_g derived with the 128³ data set, we also generated the lognormal density fields with a higher resolution of 512³ for M_s = 0.5, 1, 2.5, 5, 10, and 20. For every combination of parameters (b, M_s), we produced 30 realizations and calculated the ensemble average of σ_ln N and its statistical dispersion. Circles and error bars in Figure 4 represent the ensemble averages and the dispersions, respectively, which closely follow the curves calculated with the 128³ data set. We therefore conclude that the expressions obtained with a rather lower resolution of 128³ are also applicable to the higher resolution of 512³. The error bars obtained with the 512³ resolution indicate the ranges of variation between realizations.

In Figure 4, we overplotted the data points obtained from the HD simulations with solenoidal (blue squares) and compressive forcing (red squares) modes of Federrath et al. (2010), who provide the standard deviations of linear column density ($\sigma _{N/N_{0}}$) and logarithmic column density (σ_ln N). The filled squares show the σ_ln N values that they directly measured with the logarithmic column density and the open squares represent the values $[\ln (1+\sigma _{N/N_{0}}^{2})]^{1/2}$ calculated with the standard deviations of the column density $\sigma _{N/N_{0}}$ assuming that the column density PDF is a perfect lognormal distribution. The values estimated with the two different methods show a significant discrepancy, especially for the compressive mode, indicating that the column density PDF is not a perfect lognormal (see also Price et al. 2011 and Konstandin et al. 2012 for the discrepancies). Indeed, the PDFs for compressive forcing show large departures from a perfect lognormal function (Federrath et al. 2010). Even with this departure from the perfect lognormal function, the present results are consistent with their results within the variation ranges denoted by error bars in the figure.

In the above analyses, we have assumed HD cases. Burkhart & Lazarian (2012) fitted the MHD simulation results for solenoidal forcing with Equation (1) and obtained a best-fit A of 0.11, which is different from our result of A = 0.2. The difference might be attributable to the magnetic field effects. First, the density variance in magnetized gas is significantly lower than that in the HD counterparts (Ostriker et al. 2001; Price et al. 2011; Molina et al. 2012). The effect could be mimicked by reducing b below the value of the HD case, while assuming the same relationship between γ and M_s for the HD case. We found that the σ²_ln N–M_s relationship derived with b = 0.22 is very close to the result of Burkhart & Lazarian (2012). Second, the spectral slope of the strongly magnetized media appears to be shallower than those of unmagnetized or weakly magnetized media. Padoan et al. (2004) found that the MHD model with approximate equipartition of kinetic and magnetic energies yields a shallower density power spectrum than the super-Alfvénic model. The results of Kowal et al. (2007) also provide some indications that the power-law slopes of the supersonic, sub-Alfvénic models are slightly smaller than those of super-Alfvénic models. Their spectral slopes for the cases of a strong magnetic field are lower than the values expected from γ = 3.81 M^−0.16_s by amounts of ∼0.3–0.4. Indeed, we could reproduce the best-fit curve of Burkhart & Lazarian (2012) with a spectral slope of γ = 3.81 M^−0.16_s − 0.5 and b = 1/3. In a real situation, both effects may play roles together. Therefore, we lowered both γ and b and found that the results with γ = 3.81 M^−0.16_s − 0.3 and b = 0.29 match the result of Burkhart & Lazarian (2012), as shown by the lowest black curve in Figure 4. In the figure, the black and blue dashed lines denote the best-fit results to Equation (1) with b = 0.29 and b = 1/3, respectively.

Burkhart et al. (2010) presented γ as a function of M_s for MHD simulations with different magnetic field strengths and found that the magnetic field effect on γ is more important at lower M_s than at higher M_s. Therefore, we may need further investigations of the magnetic field effects.

Another aspect worth noting is that we can constrain the condition for Equation (1) to hold, by varying the equation between γ and $M_{{\rm {\rm s}}}$. For instance, in Figure 4, we plotted the σ_ln N–M_s relationship by assuming a constant spectral slope of γ = 2.8 (upper red curve) for b = 0.5. It is obvious that the resulting curve cannot be expressed by Equation (1). As we assume γ increases with M_s, we obtain more rapidly increasing curves than the red curve. On the other hand, if we assume that γ decreases more rapidly than γ = 3.81 M^−0.16_s with increasing M_s, the resulting curve shows a much slower increase with increasing M_s, as indicated by the lower red curve obtained for γ∝M^−0.35_s. Therefore, we conclude that the necessary condition for Equation (1) is that the spectral slope should moderately decrease as M_s increases.

We investigated the dependency of the column density variance on the 3D density variance and the power-law slope of the density power spectrum by using the fBm structures. Adopting the spectral slopes and 3D density variances that are appropriate in the HD regime, we obtained the same relationship between the column density variance and the Mach number as suggested by Burkhart & Lazarian (2012), but with a different scaling parameter A. However, when the magnetic field effects were included, the parameter A agreed with the MHD result of Burkhart & Lazarian (2012). Magnetic fields prevent the turbulent gas from reaching very low densities as well as very high densities (Molina et al. 2012), and consequently have important implications for the models of the star formation rate (e.g., Federrath & Klessen 2012).

The approach adopted in this Letter is similar to that of Padoan et al. (1997) in that they also use random density fields to derive the 3D density variance. The fBm structures also have been used to investigate the density structure of molecular clouds (Stutzki et al. 1998) and the stellar initial mass function (Elmegreen 2002; Shadmehri & Elmegreen 2012). Note that Stutzki et al. (1998) neglected the dispersion in Fourier amplitudes in generating Gaussian random fields. We compared the cases of fixed Fourier amplitudes with the results presented here, and found that the fixed amplitudes give rise to slightly stiffer spectral slopes than the present results.

The present results can be applied directly to the observations to infer the 3D density variance from the observed column density variance and spectral slope. We can vary σ_ln ρ in Equations (2) and (3) to find the best-fit value that matches the observed σ_ln N. In applying the results, we need to consider the effects of instrument noise and telescope smoothing. However, the inclusion of noise and smoothing may not severely affect the applicability of the results, as in Burkhart & Lazarian (2012).

The density fluctuation of the ISM plays a crucial role in understanding the propagation of radiation. The intensity PDF of the dust-scattered starlight in far-ultraviolet was found to be lognormal, which might be caused by the density structure of interstellar dust (Seon et al. 2011). Escape of the ionizing radiation field through the clumpy ISM has been argued to be the origin of ubiquitous diffuse ionized gas (e.g., Haffner et al. 2009). This has not, however, been clearly confirmed using realistic ISM density structures (cf., Seon 2009; Seon & Witt 2012). The present study provides a practical method to generate the realistic density structures for the radiative transfer problems.