OBSERVATIONAL DIAGNOSTICS OF SELF-GRAVITATING MHD TURBULENCE IN GIANT MOLECULAR CLOUDS

Molecular clouds are highly turbulent, magnetized and are the sites of all known star formation (Elmegreen & Scalo 2004). The details of the collapse of molecular clouds determine the key properties of the star formation rate and stellar initial mass distribution (Hennebelle & Falgarone 2012). Thus, the development of a detailed understanding of the dynamics of molecular clouds is an essential step toward a complete picture of star formation, including predicting the initial mass function.

The turbulent nature of molecular clouds is evident from a variety of observations including non-thermal broadening of molecular emission and absorption lines such as carbon monoxide (see Stutzki & Guesten 1990; Heyer & Brunt 2004) and fractal and hierarchical structures (see Vazquez-Semandeni 1994; Burkhart et al. 2013). A number of new techniques, including those studying the turbulence velocity spectrum (see Lazarian 2009, for a review) and the sonic Mach number and Alfvén Mach number (see Kowal et al. 2007; Burkhart et al. 2009, 2010, 2012; Esquivel & Lazarian 2010; Tofflemire et al. 2011), have been applied to Giant Molecular Clouds (GMCs) and also have shown that turbulence there is supersonic. In light of this, turbulence seems to play a duel role in the GMC environment of providing support on the large scales (i.e., scales of the cloud) while compressing small scales via shocks (Mac Low & Klessen 2004). Although it is clear that molecular clouds are turbulent (Elmegreen & Scalo 2004; Mac Low & Klessen 2004), magnetized (Crutcher 2012) and self-gravitating (Larson 1981; Solomon et al. 1987), the relative importance of these components is still under debate (e.g., Padoan & Nordlund 1999; Vázquez-Semadeni et al. 2008; Li et al. 2009, 2014; Kritsuk et al. 2013). As a mature theory of magnetized, supersonic, self-gravitating turbulence has yet to emerge, interpretations of observations are difficult. In light of this, numerical models have proven to be an important tool in the understanding of observations.

Several recent numerical studies have probed the observational signatures of turbulent clouds, both with and without self-gravity and with a range of Mach numbers and magnetic field strengths. Simulations, observations, and theoretic works have all pointed out the fact that compressible turbulence is important for creating filaments and regions of high density contrast (Kowal & Lazarian 2007; Burkhart et al. 2009; Federrath et al. 2010). Shocks can broaden the density/column density probability density function (PDF) and shallow the power spectral slope. The column density PDF and power spectrum (or delta variance, which is isomorphic to the power spectrum) are the chief targets for this study, due to their importance to turbulence based star formation theories (Padoan 1995; Krumholz & McKee 2005; Federrath et al. 2013). Both the density PDF and power spectrum have been shown to be sensitive to Mach number, magnetic field strength, and self gravity (Beresnyak et al. 2005; Kowal & Lazarian 2007; Kainulainen et al. 2009; Burkhart et al. 2010; Burkhart & Lazarian 2012; Collins et al. 2012; Federrath & Klessen 2013), however, less attention has been paid to higher order statistics of column density.

In this work we explore the variation of the properties of the PDFs and the power spectrum of column density as simulated clouds collapse under their own self-gravity as a function of time. The goal of this paper is to provide observers with useful methods to apply to dust extinction maps and integrated intensity maps and to explore if the PDFs/spectrum of column density behave similarly to 3D density studied in previous works such as Collins et al. (2012). In particular we are interested to explore if the collapse evolution of the cloud can be observationally discerned and pay particular attention to the time dependency of the turbulence statistics. Further, we explore two novel fits to the column density power spectrum, and use one of these techniques to disentangle the signatures of self-gravity from the turbulence and to shed light on the origin of the increase in the density power spectrum slope (Federrath & Klessen 2013). The ultimate aim of this paper is to provide set of tools to quantify the physical state of observed molecular clouds.

The layout of the paper is as follows. In Section 2 we describe the simulations studied and the codes used to generate them. In Section 3 we introduce the the column density PDF and apply it to our simulations. We investigate the evolution of the PDF power-law tail (Section 3.2), the evolution of the variance (Section 3.3), higher order moments (Section 3.4), and finally the Tsallis distribution (Section 3.5). In Section 4 we discuss the evolution of the slope of the power spectrum, a new fit to the evolving spectrum (Section 4.1); a technique to reproduce the self-gravitating spectrum from the turbulent spectrum (Section 4.2); and a technique to recover the turbulent spectrum from the self-gravitating spectrum (Section 4.3). We discuss the implications of our results in Section 5 followed by our conclusions in Section 6.

Two suites of simulations were used for this paper: one suite of non-gravitating simulations that used a fixed resolution of 512³ (which we refer to as the Godunov simulations), and one suite of self-gravitating simulations that used a 512³ root grid and four levels of adaptive mesh refinement (AMR, which we refer to as the Enzo simulations). All simulations solved the ideal magnetohydrodynamics (MHD) equations with large-scale solenoidal forcing.

The Godunov simulations used here are ten non-gravitating simulations with sonic Mach numbers ≈0.5–22 and have been used in many previous works (Cho & Lazarian 2003; Kowal & Lazarian 2007; Burkhart et al. 2009, 2010; Kowal et al. 2009, 2011). These simulations used the algorithm described in (Cho & Lazarian 2002), which uses a combination of spatially third-order essentially nonoscillatory (ENO) methods to ensure both accuracy and stability, and a third-order Runga–Kutta time integration. The Godunov models are divided into two groups corresponding to sub-Alfvénic (${B}_{\mathrm{ext}}=1.0$) and super-Alfvénic (${B}_{\mathrm{ext}}=0.1$) turbulence.

The three self-gravitating simulations were performed using the constrained transport MHD option in Enzo (MHDCT) (Collins et al. 2010; The Enzo Collaboration et al. 2013). This code uses the second order hyperbolic solver described by Li et al. (2008), the CT method of Gardiner & Stone (2005), and the divergence free interpolation method of Balsara (2001). These simulations were described in detail in Collins et al. (2012). The self-gravitating Enzo models are supersonic and super-Alfvénic, with mass chosen such that kinetic and gravitational energies are equal. They have three different values of initial magnetic field which sets the Alfvén Mach number. Throughout the text and in Table 1, we denote these different cases as low, mid, and high to stand for "low-valued," "middle-value," and "high-valued" magnetic field runs.

Table 1.  A Summary of the Simulations Presented Here

Model	${{\mathcal{M}}}_{s}$	${{\mathcal{M}}}_{A}$	Description
1	0.5	0.7	subsonic and sub-Alfvénic, Godunov
2	4.2	0.7	supersonic and sub-Alfvénic, Godunov
3	8.0	0.7	supersonic and sub-Alfvénic, Godunov
4	10	0.7	supersonic and sub-Alfvénic, Godunov
5	20	0.7	supersonic and sub-Alfvénic, Godunov
6	0.5	2.0	subsonic and super-Alfvénic, Godunov
7	4.2	2.0	supersonic and super-Alfvénic, Godunov
8	8.0	2.0	supersonic and super-Alfvénic, Godunov
9	10	2.0	supersonic and super-Alfvénic , Godunov
10	20	2.0	supersonic and super-Alfvénic, Godunov
High	8.5	6.5	supersonic and super-Alfvénic, Enzo
Mid	8.2	12.	supersonic and super-Alfvénic, Enzo
Low	9.1	30	supersonic and super-Alfvénic, Enzo

Download table as:  ASCIITypeset image

In all simulations, the ideal MHD equations are solved in a periodic box, using an isothermal equation of state ($p={c}_{s}^{2}\rho$) and a variety of sonic and Alfvénic mach numbers (${{\mathcal{M}}}_{s}=\langle | {\boldsymbol{v}}| \rangle /\langle {c}_{s}\rangle$ and ${{\mathcal{M}}}_{A}=\sqrt{\rho }\langle | {\boldsymbol{v}}| \rangle /\langle | {\boldsymbol{B}}| \rangle$, respectively). Here, ρ is density, ${\boldsymbol{v}}$ is velocity, ${\boldsymbol{B}}$ is magnetic field, p is the gas pressure, and c_s is the isothermal speed of sound. In all simulations, periodic cubes with initially uniform density and magnetic fields are driven with solenoidal forcing, using the ideal MHD equations. The self-gravitating simulations were driven until a steady state was reached, at which point gravity was turned on. Driving continued during the collapse phase. A summary of the simulations can be found in Table 1. Two representative snapshots, from the early and late times of the Enzo simulations, can be seen in Figure 1.

Zoom In Zoom Out Reset image size

For the Enzo simulations we select the Mach number, ${{\mathcal{M}}}_{s}$, virial parameter, ${\alpha }_{\mathrm{vir}}$, and mean thermal-to-magnetic pressure ratio, ${\beta }_{0}$ as

$$\begin{eqnarray}&&{{\mathcal{M}}}_{s}=\displaystyle \frac{{v}_{\mathrm{rms}}}{{c}_{{\text{}}s}}=9\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{\alpha }_{\mathrm{vir}}=\displaystyle \frac{5{v}_{\mathrm{rms}}^{2}}{3\;G{\rho }_{0}{L}_{0}^{2}}=1\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{\beta }_{0}=\displaystyle \frac{8\pi {c}_{{\text{}}s}^{2}{\rho }_{0}}{{B}_{0}^{2}}=0.2,2,20,\end{eqnarray} \tag{ 3 }$$

where ${v}_{\mathrm{rms}}$ is the rms velocity fluctuation, ${\rho }_{0}$ is the mean density, L₀ is the size of the box, and B₀ is the mean magnetic field. These can be scaled to physical clouds as

$$\begin{eqnarray}&&{t}_{\mathrm{ff}}=1.1{n}_{{\rm{H}},3}^{-1/2}\;\mathrm{Myr}\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{L}_{0}=4.6{c}_{{\text{}}s,2}{n}_{{\rm{H}},3}^{-1/2}\;\mathrm{pc}\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{v}_{\mathrm{rms}}=1.8{c}_{{\text{}}s,2}\;\mathrm{km}\;{{\rm{s}}}^{-1}\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&M=5900{c}_{{\text{}}s,2}{n}_{{\rm{H}},3}^{-1/2}{M}_{\odot }\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{B}_{0}=(13,4.4,1.3){c}_{{\text{}}s,2}{n}_{{\rm{H}},3}^{1/2}\;\mu {\rm{G}},\end{eqnarray} \tag{ 8 }$$

where ${c}_{{\text{}}s,2}=0.2\;\mathrm{km}\;{{\rm{s}}}^{-1}$ and ${n}_{{\rm{H}},3}={n}_{{\rm{H}}}/(1000\;{\mathrm{cm}}^{-3})$ are the sound speed and hydrogen number density, respectively, and we have used a mean molecular weight of 2.3 amu per particle.

In what follows we use column density maps from the Enzo and Godunov simulations to investigate the utility of often used statistical descriptors of turbulence and gravitational collapse, including the PDFs, Tsallis distribution, and Fourier power spectrum, for observations of column density, e.g., dust extinction maps. We account for fluctuations along different sight-lines relative to the orientation of the mean field in as error bars taken as the standard deviation between a measure along different sight-lines. As observers often do not have information regarding the mean field we treat the line of sight (LOS) relative to the mean magnetic field as an unknown parameter. Regarding the Enzo simulations, they are all super-Alfvénic and thus any anisotropy introduced by the mean field will be very weak (see Burkhart et al. 2014). The Godonov simulations do have sub-Alfvénic snapshots however the effect of anisotropy in these is studied in other works (see Esquivel et al. 2011; Burkhart et al. 2014).

The Probability Density Function of gas and dust in star forming regions provides important signatures of both MHD turbulence and gravitational collapse. In the case of low density MHD turbulence, where self-gravitational influences are negligible, the PDF exhibits a lognormal form for both the density and column density distributions, i.e.,

$$\begin{eqnarray}&&V(\rho )d\mathrm{ln}\rho =\displaystyle \frac{1}{\sqrt{2\pi {\sigma }^{2}}}\;\mathrm{exp}\left(\displaystyle \frac{{(\mathrm{ln}\rho -\mu )}^{2}}{2{\sigma }^{2}}\right)d\mathrm{ln}\rho ,\end{eqnarray} \tag{ 9 }$$

where $\mu =-{\sigma }^{2}/2$ is the mean of $\mathrm{ln}\rho$ and σ is the standard deviation (Blaisdell et al. 1993; Vazquez-Semadeni 1994).

Lognormal column density PDFs have been observed in various phases of the ISM including in molecular (Vazquez-Semadeni 1994; Brunt 2010; Kainulainen et al. 2011; Burkhart & Lazarian 2012; Molina et al. 2012; Kainulainen & Tan 2013) and in the diffuse warm neutral and ionized ISM (Hill et al. 2008; Burkhart et al. 2010). At the onset of gravitational collapse the shape of the turbulence induced lognormal PDF begins to become skewed toward the high density material which manifest as a power-law tail (Klessen 2000; Collins et al. 2012). Furthermore, the PDF was shown to be important for analytic models of star formation rates and initial mass functions (Padoan & Nordlund 2002; Krumholz & McKee 2005; Hennebelle & Chabrier 2008; Padoan & Nordlund 2011; Federrath & Klessen 2012, 2013).

What can be learned from the lognormal (i.e., turbulence dominated) portion of the PDF? Several authors have suggested the turbulent sonic Mach number can be estimated from the calculation of the density/column density variance (Padoan et al. 1997; Price et al. 2011) and the density/column density skewness and kurtosis, i.e., higher order moments (Kowal & Lazarian 2007; Burkhart et al. 2009).

In particular, the relationship between ${{\mathcal{M}}}_{s}$ and the variance of the logarithm of the 3D density distribution as seen in numerical models (Padoan et al. 1997; Passot and Vázquez-Semadeni 1998) generally takes the form

$$\begin{eqnarray}&&{\sigma }_{\rho /{\rho }_{0}}^{2}={b}^{2}{{\mathcal{M}}}_{s}^{2}\end{eqnarray} \tag{ 10 }$$

where ${\rho }_{0}$ is the mean value of the 3D density field, b is a constant that depends on the driving of the turbulence in question with $b=1/3$ for solenoidal forcing and b = 1 for compressive driving, and σ is the standard deviation of the density field normalized by its mean value (Nordlund & Padoan 1999; Federrath et al. 2008, 2010; Molina et al. 2012).

When taking the logarithm of the normalized density field this relationship becomes

$$\begin{eqnarray}&&{\sigma }_{s}^{2}=\mathrm{ln}(1+{b}^{2}{{\mathcal{M}}}_{s}^{2})\end{eqnarray} \tag{ 11 }$$

where $s=\mathrm{ln}(\rho /{\rho }_{0})$ and ${\sigma }_{s}$ is the standard deviation of the logarithm of density (not to be confused with ${\sigma }_{\rho /{\rho }_{0}}$).

The relationship between the observable column density standard deviation and sonic Mach number (Burkhart & Lazarian 2012) retains the same form as that of the 3D density field but with a scaling constant A:

$$\begin{eqnarray}&&{\sigma }_{\zeta }^{2}=A\mathrm{ln}(1+{b}^{2}{{\mathcal{M}}}_{s}^{2}).\end{eqnarray} \tag{ 12 }$$

The corresponding relation for the linear variance based on Equation (12) is

$$\begin{eqnarray}&&{\sigma }_{{\rm{\Sigma }}/{{\rm{\Sigma }}}_{0}}^{2}={({b}^{2}{{\mathcal{M}}}_{s}^{2}+1)}^{A}-1\end{eqnarray} \tag{ 13 }$$

where $\zeta =\mathrm{ln}({\rm{\Sigma }}/{{\rm{\Sigma }}}_{0})$ and Σ is the column density distribution with ${{\rm{\Sigma }}}_{0}$ denoting the mean value of the column density distribution.

In addition, the PDF of incremental fluctuations has been shown to be useful for studies of turbulence in the density regimes where gravity is not dominant (Esquivel & Lazarian 2010; Tofflemire et al. 2011). The Tsallis statistic provides a fit for incremental PDFs and the fit parameters describing the width and amplitude are related to the physics of the gas such as the sonic and Alfvénic Mach numbers.

In the high density regime (i.e., ${A}_{V}\gt 2$) regime, where gravity dominates the PDF shape (Schneider et al. 2013b), power-law tails form with exponent related to the star formation efficiency (Federrath & Klessen 2013). These power-law tails have been observed in numerical simulations of both density (Collins et al. 2012) and column density (Federrath & Klessen 2013) as well as observations of dust extinction in numerous star forming regions (Kainulainen et al. 2011; Schneider et al. 2013b).

In this section we focus on the column density PDFs, Tsallis function, and PDF moments in MHD turbulence simulations with and without gravity for a range of magnetic field strengths. In the case of the Enzo simulations, we pay particular attention to the time evolution in order to see how the column density PDF evolves as the cloud collapses.

3.1. Column Density PDF

Figure 2 shows the column density PDFs of the Enzo simulations for different LOS orientations (rows across) and different time steps (in units of t_ff) indicated in different colors/linestyles, with the t = 0 case indicating pure supersonic MHD turbulence plotted as a solid black line. The top row displays LOS in the x-direction (i.e., the direction of the mean magnetic field), the middle row displays the LOS in the y-direction, while the bottom row displays the LOS in the z-direction.

Zoom In Zoom Out Reset image size

We fit lognormals to the peaks of the distribution and separately fit a straight thicker red line in the power-law tail regime for $t\gt 0.35{t}_{\mathrm{ff}}$ in the column density range of ζ = 1.3–1.9. This fit range avoids the very high density portion of the PDF, which is plagued by low number statistics, as well as the peak of the distribution. We divide the PDFs into three different time regimes and point out the following visual features:

• 1.  
  At $t\leqslant 0.15{t}_{\mathrm{ff}}$, PDFs (black solid and yellow dashed lines) display general lognormal behavior. The simulation at $t=0.15{t}_{\mathrm{ff}}$ shows a wider log normal than the t = 0 simulation for roughly the same sonic Mach number. For $t\lt 0.15{t}_{\mathrm{ff}}$ in the high density range we do not observe a clear power-law tail or the tail is very steep such that it is visually indistinguishable from the lognormal.
• 2.  
  At $0.35{t}_{\mathrm{ff}}\leqslant t\leqslant 0.45{t}_{\mathrm{ff}}$, the PDFs retain the general log normal shape for log column density values $\zeta \lt 1$ regardless of LOS orientation or magnetic field strength, however for column densities $\zeta \gt 1$ power-law tails are clearly seen for time steps at $t\geqslant 0.35{t}_{\mathrm{ff}}$.
• 3.  
  At $t\gt 0.45{t}_{\mathrm{ff}}$, the PDFs (purple dotted lines) show clear power-law tails at the high column density end while the PDF distribution of low column density material has become broadened as compared with earlier times. This broadening of the low column density material has also been observed in density (see Collins et al. 2012).

Examining the PDFs of Figure 2 also shows that there are differences in the behavior for different LOS orientations and that the magnetic field strength plays a role in the PDF shape and onset of the power-law tail. In particular, turbulence with a lower magnetic field (far right column denoted with low in the title of the plot) exhibits wider PDFs at low densities and more pronounced and shallower power-law tails at the higher densities. Both the wider PDFs and shallower power-law tails are more exaggerated when the LOS is parallel (top row) to the magnetic field rather then perpendicular (middle and bottom rows). We investigate these effect quantitatively in the following subsections by exploring the PDF moments and Tsallis statistics as well as fitting the power-law tail exponential, α.

3.2. Power-law Tail

The PDF of the collapsing matter forms a power law, i.e., $V(\rho )\approx {\rho }^{\alpha }$, for densities above a critical density (Klessen 2000; Slyz et al. 2005; Vázquez-Semadeni et al. 2008). High resolution AMR simulations by Kritsuk et al. (2011) measured a range of slopes: −1.67 at intermediate densities and −1.5 at high densities. They posit that the flattening of the slope is due to the onset of rotational support, which is backed up by the analysis of the support function of that data by Schmidt et al. (2013). Recently, Federrath & Klessen (2013) investigated the power-law tail indices of self-gravitating MHD simulations and found that the high density tails are consistent with equivalent radial density profiles, $\rho \propto {r}^{-\kappa }$ with κ = 1.5–2.5. Observational constraints of the power-law tail of the column density PDF are comparable with ranges reported by simulations (Arzoumanian et al. 2011; Schneider et al. 2013b).

We plot the power-law tail exponent α versus time evolution for $t\gt 0.2{t}_{\mathrm{ff}}$ averaged over the three cardinal LOS directions with error bars indicating the standard deviation between different LOS in Figure 3. The LOS relative to the mean magnetic field is generally an unknown parameter for observational studies. We fit the slope to the range of column densities show in Figure 2 (i.e., $\mathrm{ln}{\rm{\Sigma }}/{{\rm{\Sigma }}}_{0}=1.3-1.9$). These values were chosen to avoid both the low number statistics of the high density portion of the PDF and the peak of the PDF. We will address a robust fit for the power-law tails based on maximum-likelihood fitting in a future work.

Zoom In Zoom Out Reset image size

We include the published values for the power-law tail slopes and approximate number of young stellar objects (YSOs) for five molecular clouds from Herschel data from Schneider et al. (2013a, 2014): NGC 3603, which is an active star forming regions the Auriga cloud, which is a low-mass star formation region, as well as the Orion B, Maddalena, and Aquila clouds. We discuss these clouds and their relation to the simulations in the discussion section. We omit the PDF estimates from the Carina cloud as the dynamics of Carina may include significant massive stellar feedback which is not treated in our simulations. This comparison is complementary to a study by Kainulainen et al. (2014) who showed that a proxy for the 3D density PDF could be related to the number of YSOs in several different clouds.

The power-law tail index is steepest for earlier time steps and generally becomes increasingly shallow as the cloud proceeds to collapse. This is expected from past studies of the column density power-law tail index and is also observed in the density PDF (Kritsuk et al. 2011; Collins et al. 2012; Federrath & Klessen 2012, 2013). In addition to the shallowing of the power-law index with time, Figure 3 shows that there is a strong dependency on the magnetic field. Simulations with higher field strengths show steeper slopes for all times steps due to the suppression of density enhancements by the magnetic field. The low magnetic field run shows shallower values of α across the time evolution parameter space.

The low magnetization simulation has values of α ranging from −2.5 to −0.7, which falls out of the bounds set by the observations at later times. This simulation has an Alfvénic Mach number $\approx 40$, which may be two high for realistic clouds. This may suggest super-Alfvénic clouds have Alfvénic Mach numbers in the range of 6–12, which is representative of our high and mid magnetic field simulations and has been argued for in a number of past works (Padoan & Nordlund 1999; Lunttila et al. 2008; Burkhart et al. 2009; Crutcher et al. 2009; Collins et al. 2012).

3.3. Column Density PDF Variance

We investigate the PDF variance as a function of time evolution in Figure 4 for the Enzo and Godunov simulations. We investigate several different methods of calculation of the variance in Figure 4, including directly calculating the variance of the column density distribution (${\sigma }_{{\rm{\Sigma }}/{{\rm{\Sigma }}}_{0}\mathrm{dir}}^{2},$ top panel) and directly calculating the variance of the natural logarithm of the column density distribution (${\sigma }_{\zeta \mathrm{dir}}^{2}$, center panel).

Zoom In Zoom Out Reset image size

For a distribution ${\rm{\Sigma }}={{\rm{\Sigma }}}_{i},i=1...N$ the mean value is defined as: ${{\rm{\Sigma }}}_{0}=\displaystyle \frac{1}{N}{\displaystyle \sum }_{i=1}^{N}\left({{\rm{\Sigma }}}_{i}\right)$. We define the direct calculation of the variance as

$$\begin{eqnarray}&&{\sigma }_{{\rm{\Sigma }}/{{\rm{\Sigma }}}_{0}\mathrm{dir}}^{2}=\displaystyle \frac{1}{N-1}\displaystyle \sum _{i=1}^{N}{\left(\displaystyle \frac{{{\rm{\Sigma }}}_{i}}{{{\rm{\Sigma }}}_{0}}-{\left(\displaystyle \frac{{{\rm{\Sigma }}}_{i}}{{{\rm{\Sigma }}}_{0}}\right)}_{0}\right)}^{2}.\end{eqnarray} \tag{ 14 }$$

The calculation of the variance directly from the data does not assume that the PDF follows a particular model (i.e., that it is a lognormal). We also calculate directly the variance of natural logarithm of the column density distribution as

$$\begin{eqnarray}&&{\sigma }_{\zeta \mathrm{dir}}^{2}=\displaystyle \frac{1}{N-1}\displaystyle \sum _{i=1}^{N}{\left({\zeta }_{i}-{\zeta }_{0}\right)}^{2}\end{eqnarray} \tag{ 15 }$$

Figure 4, bottom panel, shows the variance calculated by fitting a Gaussian to the peak of the distribution and measuring the variance from the fit (${\sigma }_{\zeta \mathrm{fit}}^{2}$).

The column density variance of the MHD simulations at t = 0 depends primarily on the sonic Mach number (Burkhart & Lazarian 2012) regardless of the calculation method. The distribution of gas with larger sonic Mach number shows increasing variance for both direct calculation of the variance and the variance calculated from a Gaussian fit. When gravity begins to alter the PDF from lognormal, i.e., at $t\gt 0$, the variance (${\sigma }_{\zeta \mathrm{dir}}^{2}$ and ${\sigma }_{{\rm{\Sigma }}/{{\rm{\Sigma }}}_{0}\mathrm{dir}}^{2}$) of the Enzo simulations (which have sonic Mach numbers of $\approx 10$) dramatically increases past the expectations of MHD turbulence even with Mach numbers as high as 20. This effect manifests as several orders of magnitude larger difference in the directly calculated variance (top panel) of the column density distribution. The natural logarithm of the column density distribution (middle panel) also shows a linear increase in variance as the cloud evolves with gravity. The variance as obtained from a Gaussian fit remains generally flat within the error bars across the time evolution parameter space up until about t = 0.4t_ff for the low and mid magnetic field runs and for the high field run the variance does not significantly change in the time parameter space investigated. This is because the Gaussian fit is applied to the peak of the distribution and is sensitive primarily to the lower density material dominated by turbulence and not the power-law tail high density portion of the PDF. However, there is a slight trend toward a wider lognormal at the low density end of the PDF (which can be visually seen in Figure 2). This suggests that fitting a Gaussian to the peak of the distribution can be used to dissect the turbulence dominated portions of the gas while the fraction of gas collapsing to form proto-stars may be probed with the formation of the power-law tail index (see Figure 3).

Additionally the variance of collapsing column density material shows a dependency on the magnetic field. For all three methods of variance calculation, the simulation with lower field strength shows higher values of variance as compared to the high/mid cases which are generally not distinguishable within the error bars. This effect is more pronounced and increasingly significant across different sight lines as the time evolution increases. As the magnetic field increase so does the suppression of high density clumps which causes the column density variance to be smaller compared with simulations with lower field strengths. The dependency of the column density variance on cloud magnetization in gravoturbulence simulations is somewhat surprising as the variance has been shown to be weakly dependent magnetic field strength in other works that focused only on supersonic MHD simulations without gravity (Burkhart et al. 2010; Molina et al. 2012). These results suggest the magnetic field plays a role in the global support of the cloud against gravity.

Comparison with published dust extinction column density maps from Kainulainen & Tan (2013; Table 1 of that work) shows that clouds have values of directly calculated variance ranging from ${\sigma }_{{\rm{\Sigma }}/{{\rm{\Sigma }}}_{0}}^{2}=0.25-0.64$. Taking the logarithm of these numbers for ease of comparison with the top panel of Figure 4, these values are log ${\sigma }_{{\rm{\Sigma }}/{{\rm{\Sigma }}}_{0}}^{2}=-0.6$ to −0.19, which suggests the cloud variance can not be attributed to MHD turbulence alone, and that gravity must be active up to $t=0.5{t}_{\mathrm{ff}}$ as compared with our high or mid magnetic field simulations. We over plot the Kainulainen & Tan (2013) cloud variance ranges with straight dashed lines in the top panel of Figure 4.

We compare the self-gravitating snapshots from the Enzo simulations to the column density variance–sonic Mach number relation given in Equation (13) in Figure 5. Our general conclusions indicate that two regimes exist that are defined by the importance of self-gravity:

• 1.  
  the column density variance–sonic Mach number relation without gravity: the Godunov (black squares) and t = 0 Enzo simulations (black symbols denoted in the legend) follow closely the prediction of Equation (13) given by Burkhart & Lazarian (2012).
• 2.  
  the column density variance–sonic Mach number relation with gravity: once gravity becomes important the PDF variance no longer tracks the behavior of the sonic Mach number. The more evolved the cloud's collapse, the higher the variance even for the same sonic Mach number.

Zoom In Zoom Out Reset image size

Figure 5 highlights the importance of the magnetic state of the gas in the evolution of the PDF of GMCs. The variance of the Enzo simulations depends strongly on the global strength of the magnetic field. The lower the value of the magnetic field, the higher the measured variance. This is because simulations with stronger field suppress the formation of dense cores as a function of time.

Figure 5 also compares the Enzo and Godunov column density variance–sonic Mach number relation with that of the GMCs from Kainulainen et al. (2009) as taken from Table 3 of Kainulainen & Tan (2013). The GMC variance values are all larger than the Godunov simulations and follow the upward evolution of the Enzo simulations along the variance axis. A direct comparison of the range of variance values with clouds that have ${M}_{s}\approx 10$ suggests evolution time scales of $t=0.25-0.45{t}_{\mathrm{ff}}$ over a range of magnetic field values.

3.4. Column Density PDF Higher Order Moments

In addition to the variance, the higher order moments of the PDF have also been shown in a number of works to be sensitive diagnostics of the sonic Mach number of a turbulent density field (Kowal & Lazarian 2007; Burkhart et al. 2009, 2010). However these earlier works investigated the higher order moments of the linear column density distribution (i.e., ${\rm{\Sigma }}/{{\rm{\Sigma }}}_{0}$). As shown in the previous subsections, the power-law tail of the column density distribution is a sensitive diagnostic of the evolutionary state of a collapsing cloud. A additional sensitive diagnostic that could be complimentary to the power-law tails seen in the natural logarithm of the column density distribution (ζ) are the higher order moments precisely because they characterize deviations from Gaussianity.

Skewness and kurtosis are defined by the third- and fourth-order statistical moment. Skewness is defined as

$$\begin{eqnarray}&&{S}_{\zeta }=\displaystyle \frac{1}{N}\displaystyle \sum _{i=1}^{N}{\left(\displaystyle \frac{{\zeta }_{i}-\bar{\zeta }}{{\sigma }_{\zeta }}\right)}^{3}.\end{eqnarray} \tag{ 16 }$$

If a distribution is Gaussian, the skewness is zero. Kurtosis is a measure of whether a quantity has a distribution that is peaked or flattened compared to a normal Gaussian distribution. If a data set has positive kurtosis then it will have a distinct sharp peak near the mean and have elongated tails. If a data set has negative kurtosis then it will be flat at the mean. Kurtosis is defined as

$$\begin{eqnarray}&&{K}_{\zeta }=\displaystyle \frac{1}{N}\displaystyle \sum _{i=1}^{N}{\left(\displaystyle \frac{{\zeta }_{i}-\bar{\zeta }}{{\sigma }_{\zeta }}\right)}^{4}-3.\end{eqnarray} \tag{ 17 }$$

We plot the higher order moments of the natural logarithm of the column density distribution in Figure 6. The expectations for MHD turbulence with no gravity (i.e., our t = 0 cases) are values of skewness and kurtosis bounded around zero as the lognormal distribution should have relatively small skewness and kurtosis for given sonic Mach number.⁴ An increase in sonic Mach number affects the width of the lognormal PDF and generally not the tails or peaks. Indeed, at $t\lt 0.2{t}_{\mathrm{ff}}$ the values of skewness are bounded between 0.4 and −0.4 and the values for kurtosis are bounded between −0.6 and 0.8. For $t\gt 0.2{t}_{\mathrm{ff}}$ the values of skewness and kurtosis are generally positive and increase with time evolution due to the formation of the power-law tail that creates more peaked distributions and skews the PDF tail toward the high density end. In general, we do not see a systematic sensitivity to the magnetization of the simulations, however the low and mid simulations generally show higher values of skewness and kurtosis than the simulation with high value of magnetic field.

Zoom In Zoom Out Reset image size

Figure 6 suggests that calculation of the skewness and kurtosis of the natural logarithm of the column density distribution (i.e., ζ) can be complimentary to fitting the power-law tails to the high density material. We test this idea by plotting the power-law tail index α versus the skewness and kurtosis of ζ in Figure 7. Both the higher order moment and the power-law tail index increase monotonically with time evolution and hence can be seen to increase together. Although the power-law tails have dependency on the magnetization of the gas (see Figure 3), the higher order moments of ζ are somewhat insensitive to field strength. This is clear in that the lowest magnetic field run (blue data) proceeds out to much higher values of α for the same spread in values of skewness/kurtosis. This implies that the degeneracy can be broken and that observers who measure higher order moments and power-law tails out to larger values might constrain both the energetic importance of the magnetic field as well as the time evolution of clouds.

Zoom In Zoom Out Reset image size

The many avenues of exploring PDFs, including fitting a Gaussian to the peak and a power law to the high density tail region, as well as direct calculation of the PDF moments of the entire distribution provides researchers with an avenue for separating out the contributions of MHD turbulence and gravity in the structure formation and evolution of clouds. In the next sub-section we also suggest spatial incremental PDFs to be of use toward this purpose.

3.5. Column Density Incremental PDFs: the Tsallis Statistic

In studying PDFs we study the column density field as it is given from observations. The density may be influenced by multiple processes that act differently on different scales. Thus it is advantageous to use a measure which can differentiate the properties of turbulence at different scales, e.g., at the scales of energy injection, inertial range interval, energy dissipation scale, and scales where gravitational interactions become important. To do this, it is advantageous to use incremental measures that are describe the increments of densities over the scale r, namely, $\rho (x+r)-\rho (x)$. Such incremental PDFs may be studied by different means and fitting those to the Tsallis distribution is a way to do the job.

The Tsallis distribution is often used in the context of non-extensive statistical dynamics. It was originally derived (see Tsallis 1988) from an entropy generalization to extend the traditional Boltzmann–Gibbs statistics to multi-fractal systems (such as the ISM). The Tsallis distribution has since been applied to problems over a range of applications as the distribution can describe PDFs with tails that are not exponentially bounded.

The Tsallis function of an arbitrary incremental PDF $({\rm{\Delta }}f(r))$ has the form

$$\begin{eqnarray}&&{R}_{q}=a{\left[1+(q-1)\displaystyle \frac{{\rm{\Delta }}f{(r)}^{2}}{{w}^{2}}\right]}^{-1/(q-1)}\end{eqnarray} \tag{ 18 }$$

where ${\rm{\Delta }}f(r)$ denotes the normalized PDF of incremental differences (with spatial separation r), i.e.,

$$\begin{eqnarray}&&{\rm{\Delta }}f(r)=(f({\boldsymbol{x}},{\boldsymbol{r}})-\langle f({\boldsymbol{x}},{\boldsymbol{r}}){\rangle }_{{\boldsymbol{x}}})/{\sigma }_{f},\end{eqnarray} \tag{ 19 }$$

where ${\sigma }_{f}$ the dispersion of the increments and $\langle ...{\rangle }_{{\boldsymbol{x}}}$ denotes spatial averaging over a shell of annulus of size $r=| {\boldsymbol{r}}|$. The normalization used is such that the PDF has mean value at zero and a standard deviation of unity. We denote the column density increments as: $f({\boldsymbol{x}},{\boldsymbol{r}})={\rm{\Sigma }}({\boldsymbol{x}}+{\boldsymbol{r}})-{\rm{\Sigma }}({\boldsymbol{x}})$, A normalized histogram of the incremental maps for a given lag results in our incremental PDF which is then fit with the Tsallis function.

The other parameters in Equation (18) are as follows: q is the so called "entropic index" or "non-extensivity parameter," which is related to the size of the tail of the distribution, w is a measure of the width of the PDF related to the width of the distribution, and a is the amplitude. By varying the parameter q in the Tsallis distribution it is possible to obtain distributions that range from Gaussian at $q-1$ to "peaky" distributions with large tails. The parameter q is closely related to the kurtosis (fourth order one-point moment) of the PDF, and similarly the parameter w is related to the variance of the PDF.

The Tsallis distribution reduces to the classical Boltzmann–Gibbs (Gaussian) distribution in the limit of $q\to 1$. However, for the present purpose we use it as an empiric function that fits well the properties of MHD turbulence as was shown in Esquivel & Lazarian (2010; see also Tofflemire et al. 2011). We also note that the Tsallis function was successfully used to characterize the magnetic field of the solar wind in a series of papers by Burlaga & Vinas (2004a, 2004b, 2005a, 2005b, 2006a, 2006b) and Burlaga et al. (2007a, 2007b, 2009). Esquivel & Lazarian (2010) and Tofflemire et al. (2011) used the Tsallis statistics to describe the spatial variation in PDFs of turbulent density, column density, velocity, and magnetization of MHD simulations without gravity.

Esquivel & Lazarian (2010) and Tofflemire et al. (2011) used the Tsallis statistics to describe the spatial variation in PDFs of turbulent density, column density, velocity, and magnetization of MHD simulations without gravity. Both efforts found that Tsallis provided adequate fits to their incremental PDFs and gave insight into statistics of MHD turbulence. Our present study is the first attempt to apply the Tsallis statistics to the density field obtained with the simulations that include self-gravity.

For non self-gravitating MHD simulations it was shown in Esquivel & Lazarian (2010) and Tofflemire et al. (2011) that the Tsallis fit parameters a, q, and w of the column density distribution were sensitive to both the sonic and Alfvén Mach numbers. Higher sonic Mach numbers and higher magnetic field strengths produced incremental PDFs with higher values of width (w), amplitude (a), and kurtotic index (q). In this subsection we investigate the first use of the Tsallis function to simulations of self-gravitating MHD turbulence.

Figure 8 shows the three Tsallis parameters, a, w, and q versus spatial lag (in pixels) for the Enzo simulations. Error bars are plotted by taking the standard deviation of values of the fit parameters with different LOS relative to the mean magnetic field. High, mid, and low values of magnetic field are presented for comparison in the columns going from left to right, respectively. In simulations with and without self-gravity, signs of the dissipation scales can be seen at small lag increments.

Zoom In Zoom Out Reset image size

The solid black line represents the Enzo simulation with pure supersonic MHD turbulence. Comparison of these values with Figure 8 from Tofflemire et al. (2011) show very good agreement with other supersonic MHD simulations. Once self-gravity begins to create regions of over-density, all three Tsallis parameters dramatically increase well past the expectations for supersonic turbulence.

The values of the amplitude and width of the Tsallis column density PDFs show dependencies on the magnetization of the gas. For the low and mid magnetic field simulations, the incremental PDFs are wider and have higher amplitudes compared with the high magnetization case. This is in contrast to the findings of Tofflemire et al. (2011), which found that a higher field increases the Tsallis parameters in MHD turbulence without gravity. However, when the material collapses the magnetic field acts to suppress overdense regions from forming which in turn creates incremental PDFs with lower values of width and amplitude. Within the error bars given by the LOS orientation the effect of the magnetic field is not distinguishable between the mid and low magnetic field simulations but is nearly an order of magnitude different comparing these two regimes with the high magnetic field simulation at later times (i.e., $t\gt 0.45{t}_{\mathrm{ff}}$).

The kurtotic index q shows little clear variation with changing magnetic field but does show a dependency on the time evolution of the column density. As gravity acts to create the dense clumps seen in Figure 1, the kurtotic index q is seen to increase past the turbulence only snapshot (black line). The increase in q is monotonic with increasing time step. The fact that q does not depend strongly on magnetic field (or sonic Mach number, see Tofflemire et al. 2011) but rather only on the collapse evolution suggests that this parameter might be extremely useful in breaking the degeneracy between gravity, compressibility, and magnetization in the PDF statistics of column density. We discuss this further in Section 5.

Complementary to PDFs, an essential tool for turbulence studies is the spatial density and velocity power spectrum. The turbulence energy transfer process can be studied by examining the Fourier power spectrum and the sources and sinks of energy, including the injection and dissipation scales, can be identified. The power spectrum of density and velocity (and its variants such as the structure function and delta variance) has been suggested by several authors to provide information on the spatial and kinematic scaling of turbulence, sonic Mach number, and injection/dissipation scales (Ossenkopf et al. 2001; Kowal & Lazarian 2007; Burkhart et al. 2010; Collins et al. 2012; Federrath & Klessen 2013). In this section we explore the use of the spatial power spectrum of column density maps of self-gravitating MHD turbulence for determining the dynamical evolution of clouds undergoing collapse and investigate if the power spectrum might be a complimentary tool to the PDF. We also investigate the origin of the changes in the column density power spectrum as a function of time.

The Fourier transform of the two point autocorrelation function (i.e., the spatial power spectrum) provides information on the properties of the turbulence cascade, including the injection scales and dissipation scales. The power spectrum is defined as

$$\begin{eqnarray}&&P({\boldsymbol{k}})=\displaystyle \sum _{{\boldsymbol{k}}=\mathrm{const}.}\tilde{F}({\boldsymbol{k}})\cdot {\tilde{F}}^{*}({\boldsymbol{k}})\end{eqnarray} \tag{ 20 }$$

where k is the wavenumber and $\tilde{F}({\boldsymbol{k}})$ is the Fourier transform of the field under study, which for our purposes is the synthetic column density maps.

The one-dimensional energy spectrum E(k) is related to the measured power spectrum by $E(k){dk}\propto P(k){{dk}}^{{\rm{D}}}$, where D is the dimensionality. For incompressible turbulence, the Kolmogorov power spectrum scaling (Kolmogorov 1941) in three-dimensions (3D) is $P{(k)}_{3{\rm{D}}}\propto {k}^{-11/3}$ and the energy spectrum scales as $E(k)\propto {k}^{-5/3}$. For the same $E(k)$, $P{(k)}_{2{\rm{D}}}\propto {k}^{-8/3}$, and in 1D $P{(k)}_{1{\rm{D}}}\propto {k}^{-5/3}$.

Although the Kolmogorov slope is for incompressible unmagnetized fluids, the analysis of Goldreich & Sridhar (1995) showed that the energy spectrum scaling of the incompressible cascade perpendicular to the mean magnetic field also retains the $-5/3$ slope. The Goldreich & Sridhar (1995) analysis was extended in Lazarian & Vishniac (1999) and Cho & Vishniac (2000) to include the concept of the cascade relative to the local magnetic field to obtain the correct scaling relations. Later studies confirmed the $-5/3$ slope with higher resolution simulations (Beresnyak & Lazarian 2009; Beresnyak 2012). The actual spectrum of MHD turbulence is anisotropic and scale-dependent anisotropy in the system of reference connected with the local magnetic field (see the discussion of the concept of the local magnetic field, e.g., see Cho & Lazarian 2003). The statistics of fluctuations given by the spectral slope is different when measurement are made parallel and perpendicular to the local magnetic field. However, in this paper we deal with LOS observations and in this case it is difficult to measure such anisotropy in intensity maps (see Burkhart et al. 2014 for a discussion of how this can be done in velocity centroid maps). For our purposes in this work we can use the $-5/3$ reference slope as we deal with LOS observations of super-Alfvénic turbulence and not the spectrum as measured in the local frame of reference to the magnetic field.

In the current paper we are interested in the behavior of the density/column density spectral slope in compressible self-gravitating turbulence. In the presence of supersonic turbulence, such as exists in GMCs, the density spectral slope is shallower than the relations discussed above due to shocks creating small scale enhancements of density (Beresnyak et al. 2005; Kowal & Lazarian 2007). Burkhart et al. (2010) plotted the power spectral slope of column density maps versus sonic Mach number and found that the slope of the power spectrum of ideal MHD turbulence is increasingly shallow as the Mach number increase. However they found that the power spectral slope begins to saturate toward −2 (as compared with the 2D slope of −8/3) for very high sonic Mach number, regardless of Alfvénic Mach number.

Furthermore, the inclusion of gravity in a magnetized turbulent media changes the behavior of the density spectral slope dramatically. As gravity further enhances overdense regions, the spectrum becomes increasingly shallow as more material collects on small scales (Ossenkopf et al. 2001; Collins et al. 2012; Federrath & Klessen 2013). In some cases, self-gravitating supersonic turbulence produces density structure that drive the spectral slope toward positive values. This is in contrast to non-self-gravitating turbulence where the power is dominated by large-scale structures and the power on the smaller scales is decreasing.

In this section, we investigate the evolution of the 1D column density power spectrum (which we denote as $P(k)$) as one transitions from gas dynamics dominated by supersonic MHD turbulence to self-gravitating. We also propose a model for the power spectrum of a self-gravitating fluid based on our findings. For the AMR data, the analysis was done on a coarse-grained model, where the refinement was restricted via volume-weighted average to the coarsest level. This is due to the sparse sampling of data at these higher wave numbers, which would lead to unphysical suppression of power at these wave numbers.

We plot the 1D power spectrum averaged over three lines of sight as a function of wavenumber in Figure 9 for the Enzo simulations. We overplot a solid black line in our inertial range with slope of −5/3 from log k = 1.1 to log k = 1.4 for reference. Solid colored lines are overplotted curves of Equation (21) (next subsection) with least square fit parameters given in Table 2 for each model. We find that there is no significant variation of the column density spectra with projection direction, similar to past studies such as Federrath & Klessen (2013).

Zoom In Zoom Out Reset image size

Table 2.  512³ Enzo Power Spectra Fit Parameters

Bext	t	β1	A1	β2	kc	Lc	$\langle {{\rm{\Sigma }}}_{\mathrm{max}}\rangle$	${\beta }_{{t}_{0}+\mathrm{max}(t)}$	βrec	Σc
High	0tff	−0.34 ± 0.13	0.12	−0.5 ± 0.03	N/A	N/A	4.6	−0.5 ± 0.03	−0.6	4.40826
	$0.05{t}_{\mathrm{ff}}$	−0.42 ± 0.08	0.14	−0.55 ± 0.1	N/A	N/A	4.8	−0.54 ± 0.0	−0.8	4.35521
	$0.15{t}_{\mathrm{ff}}$	−0.38 ± 0.06	0.14	−0.45 ± 0.12	N/A	N/A	9.9	−0.53 ± 0.01	−0.5	5.57656
	$0.25{t}_{\mathrm{ff}}$	−0.31 ± 0.06	0.49	−0.85 ± 0.15	120	0.241	65.7	0.39 ± 0.04	−0.4	5.41682
	$0.35{t}_{\mathrm{ff}}$	−0.36 ± 0.01	0.07	−0.12 ± 0.05	142	0.204	139.7	0.78 ± 0.02	−0.8	5.24604
	$0.45{t}_{\mathrm{ff}}$	0.1 ± 0.11	0.01	0.6 ± 0.04	700	0.041	176	0.82 ± 0.01	−0.2	5.46864
	$0.55{t}_{\mathrm{ff}}$	0.31 ± 0.07	0.01	0.78 ± 0.07	1250	0.023	261	0.86 ± 0.0134	−0.1	5.75271
	$0.60{t}_{\mathrm{ff}}$	0.41 ± 0.07	0.02	0.76 ± 0.06	1452	0.020	290	0.87 ± 0.0	−0.2	5.60185
Mid	0tff	−0.7 ± 0.05	0.12	−0.61 ± 0.03	N/A	N/A	4.8	−0.61 ± 0.03	−1.0	4.58494
	$0.05{t}_{\mathrm{ff}}$	−0.61 ± 0.11	0.09	−0.5 ± 0.05	N/A	N/A	4.6	−0.61 ± 0.0	−0.6	4.10902
	$0.15{t}_{\mathrm{ff}}$	−0.61 ± 0.04	0.16	−0.59 ± 0.08	N/A	N/A	7.4	−0.6 ± 0.0	−0.5	6.79376
	$0.25{t}_{\mathrm{ff}}$	−0.46 ± 0.07	0.62	−0.95 ± 0.03	N/A	N/A	13.6	−0.75 ± 0.03	−0.9	10.1390
	$0.35{t}_{\mathrm{ff}}$	−0.27 ± 0.0	0.38	−0.71 ± 0.04	106	0.273	80.3	0.63 ± 0.03	−1.1	5.54317
	$0.45{t}_{\mathrm{ff}}$	0.01 ± 0.16	0.03	0.26 ± 0.09	172	0.168	156	0.82 ± 0.01	−0.5	4.83240
	$0.55{t}_{\mathrm{ff}}$	0.5 ± 0.0	0.01	0.98 ± 0.17	1051	0.028	434	0.88 ± 0.0	0.3	4.60489
	$0.60{t}_{\mathrm{ff}}$	0.68 ± 0.09	0.02	0.91 ± 0.33	565	0.051	479	0.89 ± 0.0	0.1	4.05968
Low	0tff	−0.8 ± 0.01	0.21	−0.75 ± 0.14	N/A	N/A	5.4	−0.75 ± 0.14	−0.8	4.56127
	$0.05{t}_{\mathrm{ff}}$	−0.77 ± 0.1	0.17	−0.65 ± 0.06	N/A	N/A	7.1	−0.8 ± 0.0	−1.0	5.28883
	$0.15{t}_{\mathrm{ff}}$	−0.38 ± 0.04	0.15	−0.48 ± 0.12	N/A	N/A	30	−0.75 ± 0.03	−0.6	6.73403
	$0.25{t}_{\mathrm{ff}}$	−0.16 ± 0.02	0.04	0.1 ± 0.13	176	0.164	206	0.85 ± 0.01	−0.9	5.46192
	$0.35{t}_{\mathrm{ff}}$	0.78 ± 0.1	0.03	0.65 ± 0.48	350	0.083	655	0.89 ± 0.0	−0.0	4.45580
	$0.45{t}_{\mathrm{ff}}$	0.64 ± 0.1	0.05	0.79 ± 0.2	2945	0.010	784	0.89 ± 0.0	−0.7	5.62422
	$0.55{t}_{\mathrm{ff}}$	0.99 ± 0.09	0.16	0.64 ± 0.18	441	0.066	1175	0.89 ± 0.0	−0.7	12.7582
	$0.60{t}_{\mathrm{ff}}$	0.91 ± 0.03	0.24	0.79 ± 0.19	5434	0.005	1939	0.89 ± 0.0	−0.5	8.81122

Note. B_ext denotes the magnetic field regime (see Table 1); t denotes the time step in units of the free fall time; ${\beta }_{1}$ is the power spectral slope fitting from a standard power-law fit; columns 4–6 (${A}_{1},{\beta }_{2},{k}_{c}$) show the fit parameters to Equation (21); L_c shows the length scale corresponding to k_c assuming a 4 Mpc box ${L}_{c}=4\;\mathrm{pc}*2\pi /{k}_{c}$; $\langle {{\rm{\Sigma }}}_{\mathrm{max}}\rangle$ shows the maximum column density, averaged over all three axes; ${\beta }_{{t}_{0}+\mathrm{max}(t)}$ shows the slope of the turbulence with delta function, t₀ + max(t) (Section 4.2); ${\beta }_{\mathrm{rec}}$ shows the fit to the truncated power spectrum (Section 4.3); ${{\rm{\Sigma }}}_{c}$ shows the cutoff column density in code units. We can convert the code densities to physical units by multiplicative scaling factor of 1.4 × 10²² cm⁻² assuming a cloud of mean density 1000 cm⁻³ and size of 4.6 pc.

Download table as:  ASCIITypeset image

Similar to the PDFs, we find that there exists three distinct stages of evolution in the power spectrum of collapsing column density images.

• 1.  
  At $t\lt 0.15{t}_{\mathrm{ff}}$ (henceforth termed "early"), the cloud is in a purely turbulent regime, and hence the power spectrum exhibits behavior of supersonic turbulence, i.e., negative valued slopes which are shallower than the −5/3 slope.
• 2.  
  At $0.15{t}_{\mathrm{ff}}\leqslant t\leqslant 0.35{t}_{\mathrm{ff}}$ (henceforth termed "intermediate"), as the timestep increase and gravity begins to dominate the small scales (large k), the slope becomes increasingly shallow and at the largest timesteps shown, is positive in the inertial range. This is in agreement with past studies of the power spectrum of collapsing supersonic turbulence (Collins et al. 2012; Federrath & Klessen 2013). The turn over time step from negative to positive occurs around $t\approx 0.25{t}_{\mathrm{ff}}$ but is also dependent on the strength of the magnetic field.
• 3.  
  At $t\gt 0.35{t}_{\mathrm{ff}}$ (henceforth termed "advanced"), gravity dominates the power spectrum at large k values and the slopes in both the inertial rage and at larger k are positive valued.

It is interesting that the column density power spectrum (which is an observable quantity) exhibits different behavior based on time evolution. This implies that this tool maybe used on observations to not only determine the properties of turbulence in clouds but also the gravitational state of cores. We develop a functional fit to determine the turnover scale at intermediate times transitioning from turbulence dominated to gravitationally dominated in the next subsection.

We plot the power-law slopes (denoted as ${\beta }_{1}$) of the spectrum (as shown in Figure 9) versus time for the Enzo and Godunov (green symbols) simulations in Figure 10. Error bars are calculated by taking parallel and perpendicular sight lines relative to the mean magnetic field. For each sonic Mach number of the Godunov simulations, two separate Alfvénic Mach numbers are plotted.

Zoom In Zoom Out Reset image size

The Enzo simulations at t = 0 show slopes that are fully consistent with the Godunov MHD simulations for the same sonic Mach number. For $t\gt 0$ the slopes of the Enzo simulations begin to become increasingly shallow as compared with the purely turbulent scenarios, as is evident in Figure 9. However this is not significant within the error bars until $t\gt 0.2{t}_{\mathrm{ff}}$. The most significant gains in the value of the slope are made in the case of the low magnetic field. At $t\gt 0.25{t}_{\mathrm{ff}}$ this slope becomes positive and remains so as the cloud continues to evolve. For the high and mid magnetic field cases the column density power spectral slope becomes positive at $\approx t\gt 0.45{t}_{\mathrm{ff}}$. In both this work and previous works such as Federrath & Klessen (2013), positive values for the column density spectral slope are observed at evolved collapse stages.

4.1. Modeling the Power Spectrum of Self-gravitating Turbulence

Figure 9 shows that the power spectrum of the gas column density map reveals different behavior for supersonic MHD turbulence and supersonic MHD turbulence undergoing collapse. Namely that supersonic turbulence displays a negative power spectral slope while in the collapsing state the power spectrum becomes increasingly shallow past the expectations of supersonic turbulence and even becomes positive valued. This has also been confirmed in other studies, e.g., Collins et al. (2012) and Federrath & Klessen (2013).

In this subsection we address the turnover scale that occurs as the spectrum transitions from turbulence dominated to gravitational dominated in order to understand the nature of this transition and to provide observers with additional methods to disentangle the dynamics of turbulence from cloud collapse.

We propose a functional fit to the power spectrum in the form

$$\begin{eqnarray}&&P(k)={A}_{1}{k}^{{\beta }_{2}}\mathrm{exp}(-k/{k}_{c}).\end{eqnarray} \tag{ 21 }$$

This is characterized by an amplitude A₁ and power-law behavior ${k}^{{\beta }_{2}}$ which dominates the small k behavior (i.e., turbulence dominated). The scale k_c characterizes the turn over from turbulence dominated to self-gravity dominated which is observed in the power spectrum at an intermediate time evolution (see Figure 9). The advantage of the exponential form is that it can roughly describe all three stages of evolution discussed in the previous section and hence could be used by observers who do not know a priori the time scale. At t = 0 "early" timesteps the exponential can be negative to describe the dissipation of energy (a form used in previous works on MHD turbulence, e.g., Lazarian 2004). At intermediate times it can describe the turnover scale seen in Figure 9. Finally at "advanced" times Equation (21) reduces again to a power-law form as k_c becomes large. We fit the function from k = 13 to avoid the scales dominated by the turbulence driving.

Figure 9 overplots the fits of the three parameters (A₁, ${\beta }_{2}$, and k_c ) given in Equation (21) as thick solid lines. We present the fitted parameters in Table 2. For the power spectra with clear signatures of turbulence, i.e., power-law behavior in the inertial range and decreasing power in the dissipation range (large k), we only overplot the contribution of the power law in Equation (21). We note this in Table 2 under the column for k_c with N/A. Table 2 also lists the values of the slope ${\beta }_{1}$ which is calculated from a linear fit in the inertial range from k = 13 to k = 24 and plotted in Figure 10.

Figure 9 shows that Equation (21) reasonably models the turn over scale at intermediate time steps and the gravitationally dominated high k scales (i.e., small spatial scales). The fits are more robust for the high and mid magnetic field cases, whereas the low magnetic field simulation (bottom panel) is not as well fit for the mid k ranges and hence for this simulation the values of k_c are more erratic. In general, the ${\beta }_{2}$ values from the fit of Equation (21) listed in Table 2 agree within the error bar with the pure power-law fit (i.e., ${\beta }_{1}$) for early timesteps.

Table 2 shows that the turnover scale k_c generally increases with increasing timestep. This is because at $t\gt 0.5{t}_{\mathrm{ff}}$ the spectrum at the large k scales has completely transitioned from having a negative slope (turbulence dominated) to having a positive slope (gravity dominated). As k_c increases the exponential term influences Equation (21) less until again only the power-law term dominates again. At "advanced" time steps the values of k_c are equivalent to infinity since they extend beyond the range of the power spectrum being plotted (i.e., past k = 256). Thus the k_c scale of interest occurs at intermediate time steps where gravity is just starting to dominate the power spectrum. This is k_c ≈ 100–200 and occurs at timesteps from t = 0.25t_ff to t = 0.45t_ff, with the difference being attributed to the influence of the magnetic field. Using the scaling for these simulations as reported in Collins et al. (2012), this corresponds to a length scale of $L=0.005-0.2$ pc. Attention to Figure 2 shows that these timesteps are also where the power-law tails begin to form in the PDFs. Visual inspection of the column density of these snapshots reveals that this is the moment when the first cores begin to form from the turbulent cloud. Thus, the turnover scale k_c can be used as an additional diagnostic for star formation evolution in a similar manner to the PDFs and the power spectral slope.

4.2. Mimicking Self-gravity in the Turbulence Power Spectrum

The timestep of $t\approx 0.25{t}_{\mathrm{ff}}$, where gravity begin to dominate the power spectrum and PDFs of turbulent magnetized clouds, is the moment when very high density cores begin to form in our simulations. These cores exist on the smallest scales (i.e., a few pixels), yet influence the global turbulence statistics of the cloud. Thus in order to understand the affects of the small scales on the large scales a natural question arrises: can the global properties of gravity be replicated with a single (or multiple) point function which replicates a core?

We test if it is possible to replicate the power spectrum of gravoturbulence using the t = 0 snapshot of the column density (i.e., the timestep which has no gravity) and the maximum pixel value of the snapshot at some later timestep. We will call this map the "t₀ + max(t)" column density map as it is the original t = 0 map with one pixel added to it that represents the largest value of a snapshot at a later time.

Figure 11 plots the power spectrum similar to Figure 9 but overplots the power spectrum of the t₀ + max(t) column density map. In general, the correspondence is fairly good for the both the later timesteps and and earlier time steps but has trouble matching the intermediate timesteps where the turnover from turbulence dominated spectra to gravity dominated spectra takes place. At these time steps the power spectrum has a more complex behavior, i.e., can not be fit with a simple power law, which is evident by the lower value of k_c in Equation (21).

Zoom In Zoom Out Reset image size

We fit the t₀ + max(t) column density map using Equation (21) and report the slopes (denoted as ${\beta }_{\delta }$) along with the values of the maximum of the column density map in Table 2. As is evident from visual inspection of Figure 11, the largest discrepancies between the actual column density map and the t₀ + max(t) column density map are at the intermediate time steps (t = 0.25t_ff −0.35t_ff). Additionally at these time-steps the maximum value of the column density map (Table 2, column 7) increases dramatically, signaling that gravitational collapse is beginning to dominate the low k values of the power spectra of the gas. We note that our results would be similar if we had used several delta function-like points or if we had placed them in different areas of the map.

The fact that one can alter the turbulence power spectrum to mimic the gravoturbulence power spectrum provides researchers with an avenue of comparing the power spectrum of turbulence simulations with self-gravitating clouds following the prescription of maximum values given in Table 2. It also indicates that, on the scales of GMCs, the formation of cores introduce δ-function like intensity profiles to the power spectrum. This δ-function like behavior was also found in the highest density material traced by the PDFs. Finally, if it is possible to mimic the effects of the gravoturbulence power spectrum then it should be feasible to remove the signatures of gravity as well. We discuss this in the next subsection.

4.3. Restoring the Turbulence Power Spectrum in Self-gravitating Clouds

In Section 4.2 we demonstrated that it is possible to reproduce the self-gravitating power spectrum from the turbulent power spectrum, by adding a suitably chosen δ function to the turbulent spectrum. This suggests that the removal of a δ function from the self-gravitating spectrum can recover the turbulent spectrum. It was first shown by Beresnyak et al. (2005) that, for supersonic turbulence, the power the spectrum of ρ is made steeper by reducing rare density peaks, either by restriction or the logarithm. Here we try to use the relatively flat column density spectrum to recover the steeper turbulent by a suitably chosen restriction.

We define

$$\begin{eqnarray}{{\rm{\Sigma }}}_{\lt c}=\left\{\begin{array}{ll}{\rm{\Sigma }} & {\rm{\Sigma }}\lt {{\rm{\Sigma }}}_{c}\\ 0 & {\rm{\Sigma }}\gt {{\rm{\Sigma }}}_{c}\end{array}\right.\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}{{\rm{\Sigma }}}_{\gt c}=\left\{\begin{array}{ll}0 & {\rm{\Sigma }}\lt {{\rm{\Sigma }}}_{c}\\ {\rm{\Sigma }} & {\rm{\Sigma }}\gt {{\rm{\Sigma }}}_{c}\end{array}\right.\end{eqnarray} \tag{ 23 }$$

to be truncated column densities, and ${P}_{\lt c}(k),{P}_{\gt c}(k)$ to be the power spectra of each. We saw in Section 4.2 that one could construct a new field,

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{1}={{\rm{\Sigma }}}_{\mathrm{NSG}}({\boldsymbol{x}})+{{\rm{\Sigma }}}_{\mathrm{max}}(t)\delta ({\boldsymbol{x}})\end{eqnarray} \tag{ 24 }$$

where ${{\rm{\Sigma }}}_{\mathrm{NSG}}$ is the column density before the action of self gravity, and ${{\rm{\Sigma }}}_{\mathrm{max}}(t)$ is the peak column density at some time t, then the power spectrum of ${{\rm{\Sigma }}}_{1}$ is quite similar to the power spectrum of the self gravitating cloud at time t. Since it can be shown that the power spectrum of $\delta ({\boldsymbol{x}})$ is a constant with respect to k, we aim to select a cutoff column density, ${{\rm{\Sigma }}}_{c}$, that has a power spectrum that is similarly constant over a given range. Thus for some self-gravitating cloud, if ${P}_{\gt c}(k)=\mathrm{const}$ for some threshold ${{\rm{\Sigma }}}_{c}$, then we show that the spectrum of the lower column density gas, ${P}_{\lt c}$ recovers the initial turbulent spectrum reasonably well. We show this in Figure 12. In this figure, the left column is the high magnetic field simulation, and the right is the low magnetic field simulation. The top row is the full power spectra, the middle shows ${P}_{\gt c}$, which c is chosen to produce flat spectra; the bottom row shows ${P}_{\lt c}$, the truncated spectra. The black line is ${P}_{{{\rm{\Sigma }}}_{\mathrm{NSG}}}$, the initial turbulent state, and the colored lines are spectra for subsequent collapse states, as in earlier figures. To ensure ${P}_{\gt c}$ is flat we used a recursive bisection technique, varying ${{\rm{\Sigma }}}_{c}$ until the slope of ${P}_{\gt c}$ for $k\gt 10$ was near zero. While there is significant evolution in the full column density power spectrum (e.g., Figure 9), the evolution of the slope of the truncated spectrum is significantly restricted. While this technique is promising, it is is not yet able to discriminate between the initial turbulent conditions. The slope of the full spectrum, $P(k)$, for the turbulent state in the high field run, for $\mathrm{log}k$ between 1.1 and 1.4, is −0.5, while the truncated slopes vary from −0.8 to −0.2. For the mid field, the full slope is −1.0 and the variation is from −1.0 to 0.1. For the low field, the full slope is −0.8, and the truncated slope varies from −0.9 to 0.0. The last column of Table 2 shows the fits. We hope that with some refinement this technique will allow for the further separation of collapsing and turbulent clouds.

Zoom In Zoom Out Reset image size

Turbulence in GMCs is believed to be a part of a cascade that extends over 12 orders of magnitude in scale (Armstrong et al. 1995). In this paper we study the observable signatures of self-gravitating magnetized supersonic turbulence by applying the probability density functions (PDFs) and the spatial density power spectrum to synthetic column density maps generated with the Enzo code. Unlike other recent studies (e.g., Federrath & Klessen 2013) that relate density diagnostics to star formation efficiencies/rates that rely on sink particle prescriptions, we stick entirely to observable column density statistics and do not use any numerical prescriptions for star formation. In this sense we study the effects of gravitational collapse on driven turbulence and the evolution of the cloud as collapse proceeds.

We find that there exists three characterizable stages of the evolution of the PDFs and density power spectrum of the collapsing cloud which we term "early," "intermediate," and "advanced." The natural question that arises is how well do the predictions made in this paper match up with observations? The PDF of molecular clouds has been studied in a number of recent surveys and we have over-plotted measurements of the power-law tail index of those published in Schneider et al. (2013a, 2014) in Figure 3. These clouds include NGC2603, which is a high-mass active star forming regions, the Auriga cloud, which is a low-mass star formation region, and the Orion B and Aquila clouds, which are moderately star forming (Schneider et al. 2013a). Additionally, Schneider et al. (2013a) published the Herschel data for the low-mass star forming cloud called Maddalena, which they found has a power-law slope values of −3.65, however we do not include this value in Figure 3 as it is steeper than the range plotted.

How well do the ages and star formation history match with trends predicted by the simulations in Figure 3? We discuss some of the literature on each of these clouds in order of evolutionary age predicted by the PDF in Figure 3 (youngest to oldest).

• 1.  
  The Maddalena cloud has 41 young stars with disks and 33 protostars in the center of the cloud with an age estimate of a few Myr (Megeath et al. 2009).
• 2.  
  Aquila's age estimates are expected to be a few Myr (Prato & Dame 2008, p. 18).
• 3.  
  Auriga-California age estimates are not well constrained but it is suggested this cloud is not very evolved based on a high fraction of Class I and Class F YSOs (Broekhoven-Fiene et al. 2014). Harvey et al. 2013 tabulated 60 likely pre-main-sequence objects while (Lada et al. 2010) report 175-279 YSOs.
• 4.  
  Orion B has approximately 635 YSOs (Lada et al. 2010).
• 5.  
  The age estimates of NGC 3603are $10-20$ Myr (Beccari et al. 2010) making it the oldest cloud of the ones we compare with here. There are more than 10,000 stars with 7000 young stars (forming a single power-law IMF) in NGC 3603 (Harayama et al. 2008).

In light of the literature on these clouds it would seem that the PDFs are able to trace the evolutionary state of the cloud. This bodes well for future observational studies, which should combine the measures presented in this paper including the PDF variance, skewness, kurtosis, Tsallis and power spectrum to dissect the evolutionary state of clouds.

Our results also revisit the now well established debate over the impact of magnetic fields on the star formation process. The idea that magnetic fields wholly dominate the formation of stars is accepted with less enthusiasm these days. Indeed, within the dynamical picture of the interstellar medium there are factors in addition to the magnetic field that impede the gravitational collapse of a molecular cloud and make star formation less efficient, comparable to the observed rates (see McKee & Ostriker 2007). For instance, turbulent magnetic fields undergoing fast reconnection (Lazarian & Vishniac 1999), which results in reconnection diffusion (Lazarian et al. 2012), diffuse out of clouds orders of magnitude faster than the typical ambipolar diffusion timescale that is usually invoked in the traditional theory of molecular cloud evolution. Our results here show that it is wrong to disregard the role of magnetic fields completely, as we see important dependencies in the evolution of the density statistics on the magnetization of the media. Furthermore we note that Table 2 shows values of our dense core regions being within the range predicted for gravitational free fall to become faster than the reconnection diffusion timescale (i.e., around ${10}^{23}\;{\mathrm{cm}}^{-2}$, see Lazarian et al. 2012). Future works should quantify the role of reconnection diffusion in gravo-turbulent simulations such as the ones studied in this work.

Currently the gravitational state of a molecular cloud is found by comparing kinetic and gravitational energy through the virial parameter, $\alpha =5{\sigma }^{2}\;R/M$, which requires knowledge of the clouds linewidth, σ, size, R, and mass, M, as well as the implicit assumption that σ and M are spatially correlated. The techniques presented here provide a potential tool to probe the evolutionary state of a cloud using only the column density, which eliminates several sources of uncertainty in the estimation and implicitly includes information about magnetic field strength. Future numerical studies should focus on varying the sonic and Alfvén Mach numbers in order to change both the power spectral slope and dynamical importance of the magnetic field in order to determine how turbulence speeds up or impedes star formation. The filtering techniques discussed in this paper can illuminate the dynamics of turbulence even if the cloud is star forming and in an advanced state of collapse.

The current simulations aim to isolate the effects of magnetic fields and gravity on supersonic turbulence. However the use of periodic boundary conditions, solenoidal driving, and the lack of feedback are potential sources of discrepancy with real star forming clouds. Real clouds form from some sort of compressive flow, be it cloud collisions, gravitational instability, or thermal instability (Dobbs et al. 2013, and references therein) and this has been shown to impact power spectra in such simulations (Federrath & Klessen 2013). However it is possible that in the inertial range the solenoidal-to-compressive ratio reaches a universal value (Kritsuk et al. 2010). Further the periodic boundary condition imposes a conservation of total volume that may artificially enhance the amount of low density gas in the cloud. The impact of feedback is in general to inject kinetic energy, which steepens the slope by suppressing small scale structure (Sun et al. 2006). Future simulations will incorporate these effects.

Although the literature is abounding with papers regarding the 3D density PDF and its relation to the star formation efficiency, the statistical properties of velocity and magnetic field are also of vital importance to describe and understand the turbulence cascade and the star formation process. Further, these statistical properties are essential to differentiate between the "log jam" of theoretical models of star formation. Much work has gone into developing a toolbox of measurements to ascertain the physical conditions in GMCs and the ISM. The work presented here is a step toward extending this toolbox to self-gravitating turbulence. Techniques such as the Velocity Coordinate Spectrum (VCS, Lazarian 2009) can provide the injection scale of turbulence and turbulence energy density and compare the observed spectrum of clouds to analytical predictions. The studies presented here can be improved upon in the future using techniques that go beyond the PDFs and the power spectrum; for instance, earlier research with non-self-gravitating turbulence has shown that the anisotropy of velocity fluctuations can be used to find the magnetization of the medium (Lazarian et al. 2001; Esquivel & Lazarian 2005, 2010; Burkhart et al. 2014), especially if used in conjunction tools to measure the sonic Mach number (see Burkhart et al. 2013, and references therin).

It is worth summarizing the impact of the three parameters explored here (Mach number ${{\mathcal{M}}}_{s}$, Alfvén Mach number ${{\mathcal{M}}}_{{\rm{A}}}$, and collapse state) on each of the diagnostics available to us

• 1.  
  Column Density PDF power-law tails depend strongly on collapse state and ${{\mathcal{M}}}_{{\rm{A}}}$.
• 2.  
  Column Density variance from log-normal fit depends strongly on ${{\mathcal{M}}}_{s}$, slightly on ${{\mathcal{M}}}_{{\rm{A}}}$, and little on collapse state.
• 3.  
  Column Density variance from direct measurement of column density is sensitive to the collapse state and ${{\mathcal{M}}}_{s}$, and moderately correlated with ${{\mathcal{M}}}_{{\rm{A}}}$.
• 4.  
  Skewness and Kurtosis of ζ are not strongly correlated with ${{\mathcal{M}}}_{s}$; they are not strongly correlated with ${{\mathcal{M}}}_{{\rm{A}}}$, but are strongly correlated with collapse state.
• 5.  
  Tsallis fit parameters: q, w, A depend strongly on collapse state and ${{\mathcal{M}}}_{{\rm{A}}}$.
• 6.  
  Density Power Spectrum is a very strong indicator of collapse state.
• 7.  
  Velocity Power Spectrum: while not studied in this work and included here for completeness, it was shown by Collins et al. (2012) and Federrath & Klessen (2012) that the velocity power spectrum is sensitive to ${{\mathcal{M}}}_{s}$ and ${{\mathcal{M}}}_{{\rm{A}}}$, but not the collapse state. We believe that the observed low correlation of the velocity power spectrum with the collapse state is the consequence of the relatively low energy injection to the turbulence cascade during the collapse as compared to the driving energy. This matter should be explored further in the future in simulations and in observations using the VCS/VCA methods (Lazarian 2009).

The differential sensitivity of the column density statistics on the key parameters of self-gravitating turbulence opens up prospects for determining the sonic Mach number, Alfvénic Mach number and the collapse state by combining different statistics.

Turbulence, magnetic fields, and gravity are some of the key ingredients in the star formation process. Using synthetic column density maps generated with the Enzo AMR code, we investigated the robustness of the probability density functions and the spatial power spectrum for understanding and separating the roles of magnetic fields, supersonic motions and gravitational collapse in the observable column density distribution. The PDFs and power spectrum reveal three stages of cloud evolution as it progresses from diffuse turbulence dominated to collapse dominated.

Regarding the PDFs, the statistical moments and the Tsallis fit to the observable column density distribution we find that:

• 1.  
  for early times, i.e., $t\lt 0.15{t}_{\mathrm{ff}}$, the cloud has a lognormal distribution and will develop a high density power-law tail at intermediate times between t = $0.25{t}_{\mathrm{ff}}-0.35{t}_{\mathrm{ff}}$. The development of the tail, including its slope, also depends on the magnetization of the cloud. The PDF power-law tail forms earlier in clouds with lower magnetization. The tails then become increasingly shallow as the collapse proceeds with time.
• 2.  
  the directly calculated variance of the column density map (and natural logarithm of the column density map), is a sensitive diagnostic for the cloud evolution. The variance increases monotonically with time as the cloud collapses and depends on the magnetic field strength. This increase in variance with collapse causes the cloud to deviate from the sonic Mach number–variance relationship expected from non-self-gravitating turbulence.
• 3.  
  the skewness and kurtosis of the natural logarithm of the column density map are insensitive to the sonic Mach number and trace the collapse of the cloud with time, i.e., they track the formation of the power-law tail.
• 4.  
  the three Tsallis fit parameters for the incremental PDFs all strongly trace the evolution of the collapse with time: higher values correspond to more evolved cloud collapse. A strong magnetic field in a gravoturbulent cloud produces lower values of the amplitude and width of the incremental PDFs.

We find the spatial power spectrum to be complimentary to the PDFs for studies of the evolutionary state of collapsing clouds. In particular we find the following.

• 1.  
  The column density power spectrum of supersonic self-gravitating turbulence shows characteristics of a turbulence only power spectral slope at early stages of collapse ($\approx 0.15{t}_{\mathrm{ff}}$). At intermediate time-steps t = $0.25{t}_{\mathrm{ff}}-0.35{t}_{\mathrm{ff}}$, the inertial range slope becomes increasingly shallow and the dissipation range curves upwards. Eventually at advanced times the slope becomes positive and a positive sloped power law is seen down to the dissipation scales. The timescales of the changes in the slope depend on the magnetic field, with lower magnetic field facilitating earlier increases in the power spectral slope.
• 2.  
  We fit a three parameter function to the gravoturbulence power spectrum: $P(k)={A}_{1}{k}^{{\beta }_{2}}{e}^{-k/{k}_{c}}$, where A₁ describes the amplitude, ${k}^{{\beta }_{2}}$ describes the power-law behavior and the scale k_c characterizes the intermediate stage turn over from turbulence dominated to self-gravity dominated, which is observed in the power spectrum at times between $\approx 0.25{t}_{\mathrm{ff}}-0.35{t}_{\mathrm{ff}}$ depending on the strength of the magnetic field. The exponential term only becomes important (i.e., k_c becomes on the order of 100) when the first cores begin to form in our simulations. The power-law slope of the fit, ${\beta }_{2}$ is comparable to the traditional linear fit power-law slope in the inertial range.
• 3.  
  We find that the effects of self-gravity on the power spectrum can be mimicked in a turbulence only simulation by including a single point in the column density map that is equivalent to the maximum valued point in a self-gravitating map at a given time-step. This addition of a point source is reminiscent of including delta-function like behavior to the power spectrum at high k values. We provide values of the maximum points at each given time step and find that the slopes reasonably match the time-steps when k_c is very large and the exponential term to the gravoturbulence power spectra relation is negligible.
• 4.  
  We find that the effects of self-gravity on the power spectrum can be removed and the turbulence power spectrum restored through spatial filtering of the high density material.

The authors would like to thank the referee for constructive comments and suggestions which greatly improved this paper. B.B. is grateful for support from the NASA Einstein Fellowship. A.L. is supported by the NSF grant AST 1212096. B.B. and A.L. are grateful for support from the Center for Magnetic Self-organization in Laboratory and Astrophysical Plasmas (CMSO). Computer time was provided through NSF TRAC allocations TG-AST090110 and TG-MCA07S014 and the XSEDE allocation TG-AST140008. The computations were performed on Nautilus and Kraken at the National Institute for Computational Sciences (http://www.nics.tennessee.edu/) and on Stampede and Maverick at the Texas Advanced Computing Center (https://www.tacc.utexas.edu).