PDFs of molecular lines from turbulent clouds

# The Effects of Radiative Transfer on the PDFs of Molecular MHD Turbulence

Blakesley Burkhart1 , V. Ossenkopf1 ,A. Lazarian1 , J. Stutzki1 Astronomy Department, University of Wisconsin, Madison, 475 N. Charter St., WI 53711, USA Physikalisches Institut der Universität zu Köln, Zülpicher Strasse 77, 50937 Köln, Germany
1affiliationmark:
###### Abstract

We study the effects of radiative transfer on the Probability Distribution Functions (PDFs) of simulations of magnetohydrodynamic turbulence in the widely studied CO 2-1 transition. We find that the integrated intensity maps generally follow a log-normal distribution, with the cases that have best matching the PDF of the column density. We fit a 2D variance-sonic Mach number relationship to our logarithmic PDFs of the form and find that, for parameter , parameter depends on the radiative transfer environment. We also explore the variance, skewness, and kurtosis of the linear PDFs finding that higher moments reflect both higher sonic Mach number and lower optical depth. Finally, we apply the Tsallis incremental PDF function and find that the fit parameters depend on both Mach numbers, but also are sensitive to the radiative transfer parameter space, with the case best fitting the incremental PDF of the true column density. We conclude that, for PDFs of low optical depth cases, part of the gas is always sub-thermally excited so that the spread of the line intensities exceeds the spread of the underlying column densities and hence the PDFs do not reflect the true column density. Similarly, PDFs of optically thick cases are dominated by the velocity dispersion and therefore do not represent the true column density PDF. Thus, in the case of molecules like carbon monoxide, the dynamic range of intensities, structures observed and consequently, the observable PDFs, are less determined by turbulence and more-often determined by radiative transfer effects.

Radiative Transfer — ISM: structure — MHD — turbulence

## 1 Introduction

It is now accepted that the interstellar medium (ISM) of the Milky Way and other spiral galaxies is both magnetized and turbulent and that this magnetohydrodynamic (MHD) turbulence is one of the major physical processes that affect several fields of astrophysics including star formation (see McKee & Ostriker 2007 and references therein), the structure formation and evolution of different phases of the interstellar medium (Vazquez-Semadeni 1994; Ostriker et al. 2001; Elmegreen & Scalo 2004; Ballesteros-Paredes et al. 2007; Padoan et al. 2007; Chepurnov & Lazarian 2010) and magnetic reconnection (Lazarian & Vishniac 1999; Eyink, Lazarian, & Vishniac 2011). The turbulent nature of ISM111The accepted paradigm of MHD turbulence is based on the Goldreich & Sridhar (1995) model which describes the properties of Alfvénic turbulence. The decomposition in Alfvén, slow and fast modes was demonstrated in Cho & Lazarian (2003), Cho et al. (2002) and Kowal & Lazarian (2010). gases is evident from a variety of observations including electron density fluctuations, which demonstrate a Kolmogorov-type power law from kpc to AU (Armstrong et al., 1995; Chepurnov & Lazarian, 2010) and non-thermal broadening of emission and absorption lines, e.g. CO, HI (Spitzer 1979; Stutzki & Guesten 1990; Heyer & Brunt 2004).

The molecular medium, including giant molecular clouds (GMCs), represents the densest part of the turbulent cascade. Many authors have shown that molecular clouds are in fact dominated by magnetized turbulent motions, as is evident from line-width size relations, the power spectrum and the fractal nature of the observed morphology (Larson 1981; Stutzki et al. 1998). However, the exact role of turbulence in GMCs regarding the process of stellar birth is under current debate. Turbulence is a vital component in several theoretical scenarios of the star formation process such as star formation resulting from density compressions in supersonic turbulence (Elmegreen 1993; Federrath & Klessen 2012) and local core collapse in globally turbulent pressure supported clouds (Jeans 1902; Hunter 1979; Klessen et al. 2000). A recently proposed model of magnetic flux transport out of clouds, namely, the model of “turbulent reconnection diffusion” (Lazarian 2005; Lazarian et al. 2010; Santos-Lima et al. 2010, 2011) is based on the theory of fast reconnection of magnetic fields in a turbulent medium (Lazarian & Vishniac 1999).

Despite the recognition that turbulence (particularly supersonic turbulence) and magnetic fields are important to studies of molecular clouds and star formation, the community still faces substantial challenges in studying and characterizing the MHD properties in a quantitative way. This is, in part, because the fundamental parameters of interstellar turbulence, such as the sonic and Alfvénic Mach numbers, as well as the virial parameter and turbulence driving mechanisms (i.e. the forcing parameters), are difficult to determine due to the challenges of imaging small scale turbulence (Bertoldi & McKee 1992; Federrath et al. 2010; Gaensler et al. 2011; Burkhart et al. 2012a) and because measurements of magnetic fields are are still marginal and time consuming.

In addition to observational data, the analysis of simulations have proven to be very useful for studying turbulence in the diffuse ISM and in star forming regions. The most common statistical tool used for studies of turbulence is the spatial power spectrum of column density maps or electron scintillation and scattering measurements (Armstrong et al. 1995; Elmegreen & Scalo 2004) although more recently several studies have undertaken measurements of the velocity power spectrum (Lazarian 2009). More recently, the power spectrum of column density maps has also been used to estimate the star formation efficiency rates in numerical simulations (Federrath & Klessen 2013). Although the power spectrum is useful for obtaining information about energy transfer over scales, it does not provide a full picture of turbulence, partially because it only contains information on Fourier amplitudes. In light of this, many other techniques have been developed to study and parametrize observational magnetic turbulence. These include higher order spectra, such as the bispectrum, higher order statistical moments, topological techniques (such as genus), clump and hierarchical structure algorithms (such as dendrograms), principle component analysis, Tsallis function studies for ISM turbulence, and structure/correlation functions as tests of intermittency and anisotropy (for examples of such studies see Heyer & Schloerb 1997; Chepurnov et al. 2008; Burkhart et al. 2009; Chepurnov & Lazarian 2009; Kowal, Lazarian, & Beresnyak 2007; Goodman et al. 2009a; Esquivel & Lazarian 2010; Roman-Duval et al. 2011; Tofflemire, Burkhart, & Lazarian 2011; Burkhart et al. 2012b). Wavelets methods, and variations on them such as the -variance method, have also been shown to be very useful in characterizing the distribution of inhomogeneities in data (see Ossenkopf & McLow 2002; Ossenkopf, Krips, & Stutzki 2008a, 2008b; Schneider et al. 2011).

However, many of the statistical studies being conducted fail to take into account key radiative transfer properties and observational constraints, such as the varying optical depth, when investigating the structure and statistics of turbulence. This is particularly true for star forming molecular clouds, as absorption effects are important for all abundant molecular tracers, like e.g. CO and CO. Furthermore, recent studies have brought into question whether the emission of CO or its isotopes can provide a good representation of the true density distribution in molecular clouds and under which conditions do they really trace H (see Goodman et al. 2009b; Glover et al. 2010; Shetty et al. 2011a,b). Four different effects may change the so-called -factor, i.e. the ratio between the CO intensity and the column density of molecular () gas. These are variations of the abundance of CO relative to H, the temperature of the molecular cloud, the velocity range, and density effects: subcritical collisional excitation in thin gas, and saturation due to high optical thickness for large column densities. Previous systematic studies (Shetty et al., 2011a, b) combined all effects but focused on the impact of a changing molecular abundance. To obtain a more complete description in terms of the impact of the different physical effects on the measured structure of a molecular cloud, we investigate here in detail the impact of sub-thermal excitation and radiative transfer effects on the emission of turbulent molecular clouds. We use isothermal MHD simulations to simulate molecular line maps of the CO 2-1 transition from the models comparable to the observed maps. We generate the molecular line maps using the Ossenkopf (2002) radiative transfer algorithm (SimLine-3D). We then investigate the resulting CO line mps, in particular their Probability Distribution Functions (PDFs) and their descriptors, such as variance, skewness, and kurtosis and the Tsallis function fit, and study how these quantities change when changing the optical depth, density, and abundances (see Table 1 for details on our parameter space).

For this paper, we concern ourselves only with the 2D quantities (column density, integrated intensity) that are available from direct observations. We will denote the 2D column density or integrated intensity distribution as . We will use to stand for the logarithm of the 2D quantities, (i.e. , where is the mean value of ).

The paper is organized as follows. In § 2 we describe the PDF methods used in this paper. In § 3 we describe the simulations and outline the main points of the radiative transfer algorithm. In § 4 we discuss the use of PDFs on turbulent clouds and explore the PDFs of our models. In § 5 we discuss the variance vs. sonic Mach number relationship for our models. In § 6 we discuss the relationship between the skewness and kurtosis and the sonic Mach number for our models. In § 7 we apply the incremental PDF functional fit known as the Tsallis distribution to our models and investigate the relationship the fit parameters have with the optical depth, sonic and Alfvénic Mach numbers. Finally, in § 8 we discuss our results followed by the conclusions in § 9.

## 2 Probability Distribution Functions as a Tool for Turbulence Studies

PDFs are one of the statistical tools developed to study turbulence in both diffuse and self-gravitating clouds. PDFs and their quantitative descriptors have been used to study turbulence in a variety of astrophysical contexts including the diffuse ISM (Hill et al. 2008; Berkhuijsen & Fletcher 2008), the solar wind (Burlaga et al. 2004a) and the molecular ISM (Padoan et al. 1997; Kainulainen & Tan 2012).

Several authors have suggested that PDFs of column density data can be used to predict important parameters of ISM turbulence in the diffuse and molecular medium such as the sonic Mach number (Burkhart et al. 2009, 2010; Price, Federrath, & Brunt 2011), the star formation rate in GMCs (Krumholz & McKee 2005; Federrath & Klessen 2012), and the initial mass function (Hennebelle & Chabrier 2008, 2011). Two predominant ways of exploring the PDFs of turbulence have been investigated: one focuses on the statistics of linear column density while the other takes advantage of the generally log-normal nature of turbulent fields. In either case, researchers have either calculated the moments of the PDF distribution such as the variance, skewness, and/or kurtosis or fitted a function to the PDF and used the resulting fit parameters as descriptors of turbulence (e.g. the Tsallis function, Esquivel & Lazarian 2010; Tofflemire, Burkhart, & Lazarian 2011).

In regards to the moments, several studies have numerically investigated the relationship between the 3D density PDF width (i.e. the variance or standard deviation) and the root mean square (rms) sonic Mach number for linear and natural-log density (Padoan et al. 1997; Price, Federrath, & Brunt 2011; Molina et al. 2012). For a distribution the mean value and variance are defined as: and . The primary advantage of taking the logarithm of the density field is that turbulent density fields are generally log-normal, and thus researchers can take advantage of analytically relationships of log-normal distributions. A key assumption in these models is that there is a linear relationship between the variance in linear density and the sonic Mach number, i.e,

 σ2ρ/ρ0=b2M2s (1)

where is the mean value of the 3D density field, is a constant of order unity and is the standard deviation of the density field normalized by its mean value. When taking the logarithm of the normalized density field this relationship becomes:

 σ2s=ln(1+b2M2s) (2)

where and is the standard deviation of the logarithm of density (not to be confused with ). Reported values for vary widely for different numerical studies and under different turbulent driving schemes. Values of from numerical studies are found to be between 0.26 and unity. Generally, the value depends on the driving of the turbulence in question with for solenodial forcing and for compressive driving (Nordlund & Padoan 1999; Federrath, Klessen, & Schmidt 2008).

A key limitation of many of the variance-Mach number studies mentioned above is that the relationships are derived only for 3D density, which is not observable. Recently, several authors have attempted to solve this problem. A study by Brunt, Federrath, & Price (2010) devised a method to reconstruct the 3D variance density from the projected 2D variance, however the study found a error for different sight-lines for the case of ideal MHD turbulence and did not include effects of self-gravity or radiative transfer (which can significantly alter the distribution of integrated intensity). A study by Burkhart & Lazarian (2012) modified Equation 2 to be directly applied to 2D column density maps.

The variance is not the only PDF descriptor that can be related to the sonic Mach number. Studies of the higher order moments of density and column density (for linear distributions) by Kowel et al. (2007) and Burkhart et al. (2009, 2010) showed that the skewness and kurtosis have a nearly linear relationship with the sonic Mach number. The 3rd order moment or skewness is defined as:

 γX=1NN∑i=1(Xi−¯¯¯¯¯XσX)3 (3)

and the 4th order moment kurtosis as:

 βX=1NN∑i=1(Xi−¯¯¯¯¯XσX)4−3 (4)

Skewness describes the asymmetry of a distribution about its mean. The fourth order moment, kurtosis, is a measure of a distribution’s sharpness or flatness compared to a Gaussian distribution. Skewness and kurtosis are therefore unitless numbers which require no normalization between simulations/observations (variance, however, is scale dependent). They also require no assumption of log-normality for the PDF in question. As the sonic Mach number increases, shocks create density enhancements and hence the skewness and kurtosis increase. This reflects a departure from Gaussian subsonic distributions, which have skewness/kurtosis 0.

Additional PDF descriptors used for studies of ISM turbulence include the Tsallis function. The Tsallis distribution is a function that can be fit to incremental PDFs of turbulent density, magnetic field, and velocity. One of the Tsallis fit’s first practical application for studies of turbulence were on PDFs of the solar wind. Burlaga and collaborators applied the Tsallis distribution to the solar wind magnetic and velocity field data (Burlaga & Viñas 2004a, 2004b, 2005a, 2005b, 2006) to describe the temporal variation in PDFs measured by the Voyager 1 & 2 spacecrafts. Esquivel & Lazarian (2010) and Tofflemire, Burkhart, & Lazarian (2011) used the Tsallis statistics to describe the spatial variation in PDFs of density, velocity, and magnetic field of MHD simulations similar to those used here. Both efforts found that the Tsallis distribution provided adequate fits to their PDFs and gave insight into the sonic and Alfvénic Mach numbers of MHD turbulence.

To our knowledge, no numerical study has been performed to address the effects of radiative transfer and optical depth on the above mentioned PDF descriptors and their relations to the sonic Mach number.

## 3 MHD Simulations and Radiative Transfer

We generate 3D numerical simulations of isothermal compressible (MHD) turbulence by using the the Cho & Lazarian (2003) MHD code and varying the input values for the sonic and Alfvénic Mach number.

We briefly outline the major points of their numerical setup (for more details see Cho et al. 2002, Cho & Lazarian 2003).

The code is a second-order-accurate hybrid essentially nonoscillatory (ENO) scheme (Cho & Lazarian 2003) which solves the ideal MHD equations in a periodic box:

 ∂ρ∂t+∇⋅(ρv)=0, (5) ∂ρv∂t+∇⋅[ρvv+(p+B28π)I−14πBB]=f, (6) ∂B∂t−∇×(v×B)=0, (7)

with zero-divergence condition , and an isothermal equation of state , where is the gas pressure. On the right-hand side, the source term is a random large-scale driving force (). We drive turbulence solenoidally, at wave scale equal to about 2.5 (2.5 times smaller than the size of the box). This scale defines the injection scale in our models in Fourier space to minimize the influence of the forcing on the generation of density structures. Density fluctuations are then generated by the interaction of MHD waves. The time is in units of the large eddy turnover time () and the length in units of , the scale of energy injection. The magnetic field consists of the uniform background field and a fluctuating field: . Initially .

We divided our models into two groups corresponding to sub-Alfvénic () and super-Alfvénic () turbulence. For each group we compute several models with different values of gas pressure, which is our control parameter that sets the sonic Mach number (see Table1, second column). We run compressible MHD turbulent models, with 256 and 512 resolution, for crossing times, to guarantee full development of the energy cascade. Since the saturation level is similar for all models and we solve the isothermal MHD equations, the sonic Mach number is fully determined by the value of the isothermal sound speed, which is our control parameter. The models are listed and described in Table 1.

For our tests we choose a total cube size of 5pc and a gas temperature of 10K. The cube is observed at a distance of 450pc with a beam FWHM of 18” and a velocity resolution of . We note that the original simulations set the sound speed by varying the isothermal pressure and holding the RMS velocity constant. However, now we set a constant temperature of 10K for all the simulations. In this case, it makes sense to hold the sound speed constant and vary the RMS velocity in order to set the isothermal sonic Mach number. However, as we show in the appendix, the resulting PDFs are unaffected by the rescaling of the velocity field, and the parameters of critical importance remain the optical depth and sonic Mach number. Additionally, the radiative transfer has a dependency on the temperature, however here we concern ourselves with the density regime of GMCs that is known to be approximately isothermal, i.e. at intermediate densities, such as the norm case in Table 1, which are larger than the transition density to the CNM but smaller than the densities leading to star formation.

We vary both the number density scaling factor (in units of , denoted with the symbol ) and the molecular abundance (CO/H, denoted with the symbol ) in order to change excitation and optical depth . To represent the typical parameters of a molecular cloud we choose standard values for density, 275 , and abundance, (this is what we will refer to as the normalized or norm parameter setup), providing a slightly optically thick CO 2-1 line at the highest densities with line center optical depths of a few ((norm) in Table 1 in column 5). To investigate the impact of line saturation and subthermal excitation we vary density () and abundance () by factors of 30 to larger and smaller values compared to the (norm) values. We list all the parameters of the complete model set in Table 1.

We show examples of our synthetic CO 2-1 integrated intensity maps for a range of optical depths, abundances (), and initial densities () for a supersonic case (model 6) and a subsonic case (model 2) in Figure 1 . Both cases are sub-Alfvénic.

Increasing density or abundance (far right images in panels A and B of Figure 1) leads to similar effects of optically very thick lines with a peak optical depth around . The map is dominated by the velocity dispersion, not the true column density structure. For lowered abundances (bottom left images in panels A and B of Figure 1), we get a more accurate representation of the true column density structure. For lowered densities (top left images in panels A and B of Figure 1), the regions of low density become invisible in CO because the molecule is no longer thermally excited.

Examination of Figure 1 clearly show that different sonic Mach numbers provide different distributions of the integrated intensity maps for CO 2-1 emission. However, comparison between different ranges of average density/abundance ratio (and as a result, different values of optical depth) also show differences in the distribution and intensity values of the integrated intensity maps (and therefore the resulting column density maps of for a given X factor). In fact, in many cases presented in Figure 1 (especially the highly optically thick cases), it is clear that the radiative transfer dominates over turbulence in terms of the distribution and structures present. Thus for carbon monoxide, the common statistics of turbulence studies must be revisited to include these effects. In the next sections we will investigate the effects of optical depth on the PDF moments and incremental PDF Tsallis fits in order to determine over what ranges of optical depth these statistics are still sensitive to the Mach numbers for CO 2-1 emission.

## 4 PDFs of 13CO Turbulent Maps

We calculate the PDFs for all models listed in Table 1 for both the linear distribution and taking the natural logarithm of the maps. We make plots of the PDFs of the models for super-Alfvénic turbulence. We plot the linear PDFs (i.e. PDFs of ) of the column density (red lines) and integrated intensity in Figure 2. The panels are ordered subsonic to supersonic going left to right. Different colored lines and symbols represent our varying radiative transfer parameter space. For the subsonic case (far left panel) it is clear that none of the models that include radiative transfer exactly match with the subsonic column density PDF of , which is due to the fact that mutual radiative pumping takes over the role of thermalization of the lines. As the sonic Mach number increases (center panel with and far right panel with ) we can also see a general departure from the column density PDF as the radiative transfer environment alters the integrated intensities. Additionally, it is clear that the ’norm’ (yellow asterisk) model PDF with fits very closely to the column density PDF.

The column density PDF becomes increasingly skewed and kurtotic as the sonic Mach number increases compared to the Gaussian-like distribution seen in the subsonic panel. This is similar to findings from past works such as Kowal et al. (2007) and Burkhart et al. (2009, 2010), who suggested that these trends are due to high density regions created by shocks which skew and peak the distribution away towards higher density values. However, neither of these studies explored distributions of integrated intensity maps. As the optical depth decreases, either due to decreasing the abundance (purple squares) or the average density (blue diamonds), the distribution becomes increasingly skewed and kurtotic compared with the column density. This is because the emission is concentrated only in the regions of higher density and abundance while the thin gas is not strongly emitting, creating low emission holes. For cases of higher optical depth, either from increasing the density (black plus signs) or increasing the abundance (green triangles) the distributions are significantly less skewed and peaked compared with the column density. This is due to the effects of self-absorption in the gas. The cloud becomes optically thick before it can be fully sampled, thus masking out any shock enhancements that could be seen in the cloud. Thus in the case of high optical depths, the distributions are similar to the low Mach-number cases.

We make plots of the PDFs for the logarithm of the distribution () of the models for super-Alfvénic turbulence in the top plot of Figure 3. We plot the PDFs of the logarithm of column density (red lines) and logarithm of integrated intensity with the same models as in Figure 2.

All distributions appear approximately log-normal, however some deviations can be seen in the subsonic PDFs (far left panel). With line saturation, the PDF can be expected to be asymmetric, with a tendency to pile up at larger intensities. This effect was also observed in Shetty et al. 2010. For the subsonic case, all models with radiative transfer show slightly larger widths than the true column density, regardless of optical depth.

In these supersonic cases (which all generally appear log-normal), the norm case with a mixed optically thin/thick medium (orange stars) has a width that is generally the closest to the column density. The two optically thin cases (blue diamonds and purple squares) have very large widths and begin to deviate from log-normal. This is again due to quenching of the emission due to a limited range of densities/ abundances, which in turn increases the variance in the distribution due to larger contrast between the few areas with enough density/tracer material and the invisible gas. The optically thick cases (green triangles and orange stars) have smaller PDF widths than the true column density.

Figures 2 and 3 show that the PDFs of the 2D distribution of MHD turbulence as seen from observations will not only depend on the sonic Mach number, but also on the radiative transfer environment and the optical depth of the medium. This fact has not been accounted for in previous studies, and may affect the derived relationships between the sonic Mach number and the variance, skewness, and kurtosis of linear and logarithmic column density/integrated intensity. It may also explain the discrepancy and uncertainties between the PDFs of numerical simulations to those of observations (see Goodman et al. 2009b; Price, Federrath & Brunt 2011; Burkhart et al. 2009,2010; Brunt 2010). We address these issues in the next two sections.

## 5 Ms-σ Relation For Column Density and Integrated Intensity Maps

In this section we investigate the relationship between the root-mean-squared (RMS) sonic Mach number and the variance of both the linear column density/integrated intensity and logarithmic column density/integrated intensity. The relationship between the 3D variance and sonic Mach number is given by Equation 1 for density, with parameter being empirically derived from numerical simulations.

More recently, Burkhart & Lazarian (2012) formulated a variant of Equations 1 and 2 for 2D turbulent maps as:

 σ2Σ/Σ0=(b2M2s+1)A−1 (8)

and

 σ2ζ=A×ln(1+b2M2s) (9)

,where (see Federrath et al. 2010) and were found from a best fit. However, they did not investigate the effects of radiative transfer that are critical for application to the molecular medium. In the following subsections, we will investigate the robustness of Equations 8 and 9 on our CO integrated intensity maps.

### 5.1 Linear Variance vs. Sonic Mach Number

In Figure 4 we plot the linear variance of () vs. sonic Mach number for all our simulations. The dashed lines show the the best fit given by Equation 8 including radiative transfer effects and using b=1/3 for solenodially driven simulations (see Federrath et al. 2010). Error bars are created by taking the standard deviation between the variance of integrated intensity maps with line-of-sight parallel and perpendicular to the mean magnetic field. The resulting parameter is given in Table 2.

Inspection of Figure 4 reveals that values of , for sonic Mach numbers greater than 2, are clearly dependent on the radiative transfer environment. When we include the effects of radiative transfer, the closest match to the column density is surprisingly not for optically thin tracers, but rather our norm case with . Equation 8 provides a good fit to the true column density (red symbols), the mixed optically thin/thick case (yellow line) and the high optical depth cases (green and black lines) at supersonic Mach numbers. For subsonic Mach numbers, the global radiative coupling for Mach numbers around unity and below distorts the relationship between the sonic Mach number and variance.

### 5.2 Logarithmic Variance vs. sonic Mach number

We analyze the variance of the logarithmic integrated intensity maps and their relationship with the sonic Mach number (i.e. we plot vs. ) in Figure 5. Error bars are created by taking the standard deviation between the variance of integrated intensity maps with line-of-sight parallel and perpendicular to the mean magnetic field. The dashed lines represent Equation 9 when using and the parameters from the fit of the linear distributions in Table 2. The resulting fits are very similar to what was found for the case of the linear variance vs. sonic Mach number. Equation 9 provides an acceptable fit to the true column density (red symbols), the mixed optically thin/thick case (yellow line) and the high optical depth cases (green and black lines) at supersonic Mach numbers. For subsonic Mach numbers, again we find that the global radiative coupling distorts the relationship between the sonic Mach number and variance. In general, values for are larger in cases where the optical depth is low and smaller for higher optical depth.

## 6 Ms-Skewness and Kurtosis of Column Density and Integrated Intensity Maps

The higher order moments skewness and kurtosis (definitions are given in Equations 3 and 4) have also been shown to have a strong dependency on the sonic Mach number of both 3D density and 2D column density (Kowal et al. 2007; Burkhart et al. 2009, 2010).

We calculate the higher order moments for the linear distribution of the integrated intensity maps and for the true column density as per Equations 3 and 4. We plot the skewness and kurtosis vs. sonic Mach number for the integrated intensity maps in Figures 6 and 7, respectively. For the kurtosis vs. sonic Mach number, we over-plot a black dashed line to show the linear fit used in Burkhart et al. 2010).

For both Figures 6 and 7, in the case of optically thin media (blue and purple), the higher order moments increase more dramatically (i.e. non-Gaussianity increases) with increasing sonic Mach number than the true column density. In particular, the low density case (blue diamonds, model ) departs from the general monotonic trend seen in the true column density and other models. The low density case also has extremely large values of kurtosis (greater than 60) for our highly supersonic models, and we do not show them on the plot in order to preserve the visualization for the other models. The model that most closely matches the true column density in both skewness and kurtosis is the norm case with . As the optical depth increases, the dependency of the sonic Mach number on the skewness and kurtosis becomes less apparent. In both cases of changing the density and abundances it is clear that the higher order moments are very sensitive to the radiative transfer environment in addition to the sonic Mach number.

We chose here to focus on the relationship between the moments of the column density distribution and the sonic Mach number for the CO integrated intensity maps. In the above analysis, it is clear that an important parameter to take into account when applying PDF-Mach number relations on CO is the optical depth. In this case we also show how the moments change with varying in Figure 8 for our sub-Alfvénic models. If the sonic Mach number is known (e.g. via velocity dispersions), then one can potentially estimate the order of magnitude of the optical depth in supersonic clouds.

## 7 Tsallis Statistics

A separate avenue for investigation of PDFs of interstellar turbulence has been to use incremental PDFs along with the Tsallis function. The Tsallis distribution (Tsallis 1988) was originally derived as a means to extend traditional Boltzmann-Gibbs mechanics to fractal and multi-fractal systems.

The Tsallis function is used to fit incremental PDFs, i.e. PDFs of , where refers to a spatial average. We note that the increments do not need to be spatial but could also be temporal increments.

The Tsallis function of an arbitrary incremental PDF has the form:

 Rq=a[1+(q−1)Δf(r)2w2]−1/(q−1) (10)

The fit is described by the three parameters , , and . The parameter describes the amplitude while is related to the width or dispersion of the distribution. Parameter , referred to as the “non-extensivity parameter” or “entropic index”, describes the sharpness and tail size of the distribution.

The Tsallis fit parameters are in many ways similar to statistical moments explored in the previous sub-section. In regards to the Tsallis fitting parameters, the parameter is similar to the second order moment variance while is closely analogous to fourth order moment kurtosis. Unlike higher order moments, however, the Tsallis fitting parameters are dependent least-squares fit coefficients and are more sensitive to subtle changes in the PDF.

We apply the Tsallis fit to PDFs of increments using the CO 2-1 integrated intensity maps varying the average density and abundance factor. We show an example of the Tsallis fits for super-Alfvénic CO 2-1 turbulence in Figure 9. We find that the Tsallis function fits incremental PDFs of integrated intensity maps well, including both regimes of optically thin and thick. As an example, we only show the results of varying the abundance in Figure 9 and omit varying the density parameter space. The blue line is the optically thin case, the red line represents our standard case with and and the black line is the most optically thick case with abundance value larger than the norm. Columns show different sonic Mach number while rows show different incremental lag values. Solid lines represent the Tsallis fit while the dotted lines are the PDFs.

We find that the incremental PDFs are less peaked for cases of high lag, high optical depth and low sonic Mach number. Examination of Figure 9 shows that, for the fit values of the width and amplitude of the PDF of the =7.0 case shows a strong dependency on the optical depth/ radiative transfer environment over all lag values. Interestingly, we see less of a dependency on the radiative transfer for the =2.0 and =0.4 cases.

We further investigate the dependencies of the incremental PDF width and amplitude on sonic Mach number and optical depth by plotting the log of the width () and log of the amplitude parameter () vs. spatial lag for cases of varying density and abundance in Figure 10. Rows show varying sonic Mach number while sub- and super-Alfvénic turbulence is shown left and right. Cases with higher magnetization (lower Alfvénic number) generally have higher values of the width and amplitude. The symbol and color scheme used is similar to all previous plots used in this paper.

The log of the amplitude parameter is shown in the top panel of Figure 10. The =7.0 case shows the clearest distinction between cases of high and low optical depth over all lags. The amplitude is highest for the cases of high sonic Mach, especially when the medium is optically thin. Once the optical depth has increased into the optically thick regime, differences between sonic Mach number become hard to distinguish. When the medium is supersonic and the average density is low (blue line, top row) the amplitude of the incremental PDF is much larger than the high abundance case, which implies the Tsallis function is also sensitive to changes in density/abundance rather than just the optical depth.

The width parameter (in the bottom panel of Figure 10) shows similar trends as seen in the amplitude parameter. The =7.0 case shows the clearest distinction between the varying radiative transfer parameter space. For supersonic turbulence, the widths are highest in the presence of an optically thin medium. Once the optical depth has increased into the optically thick regime, differences between sonic Mach number become hard to distinguish, however, there still remain some differences in Alfvénic Mach number.

For both the amplitude and width parameter, the case that most clearly follows the true column density is the norm case with . For this case, differences in both sonic and Alfvénic Mach number can be seen. This makes the Tsallis function particularly useful for clouds whose optical depth are around unity. Furthermore, if the sonic Mach number is known and is found to be supersonic, the Tsallis fit can be used to determine the radiative transfer parameter space.

## 8 Discussion

This study stresses the importance of radiative transfer effects in studying the properties of the observable line intensity distribution of the turbulent ISM. They strongly affect the molecular line maps in turbulent clouds. We restricted the analysis to the PDFs of the CO J2-1 transition. We expect the same qualitative behavior with higher transitions rather following the reduced density case and CO following the high abundance case. While we focus our study on carbon monoxide, our conclusions may well be applicable to other species, e.g HI.

This study is complementary to a study done by Shetty et al. (2011a,b), which found that the X factor depends strongly on the average density and CO abundances. As was observed with the PDFs in this work, high values of volume density and/or CO abundance can cause the CO line to saturate, which creates integrated intensity maps that miss important structure in the true column density map. In the case of optically thin CO and low density, the lowest CO density/abundance regions will not collisionally excite, and again we miss out on important structure in the true column density map. This harkens back to a study by Goodman et al. (2009b), who showed that one of the biggest downfalls of carbon monoxide tracers is its limited dynamical range. Goodman et al. (2009b) pointed out that dust, while more limited in terms of resolution, has a much larger dynamical range with which to trace out the true column density map. However, it cannot trace any velocity information through line profiles.

In this study, we find that CO is able to trace the PDF statistics of the true column density best for case with using average interstellar values for the density () and abundances (, Burgh et al. 2007), as represented in our norm case.

It may seem counter-intuitive that the best tracer of the true column density PDF is the one with unity, since many past studies have always opted for the most optically thin observations available in order to avoid saturation effects. However, it seems the key issue is the limited dynamical range of CO. Figure 1 shows that the dynamical range is higher in the cases of the optically thick cases (due to the fact that these cases are dominated by dispersion). However, the case with optical depth on order unity (central panels of Figure 1 shows that this case still traces most of the true column density structure while having the largest dynamic range. In this case, getting the true column density can no longer be considered only in terms of obtaining the correct X factor (which has its own dependencies on abundance, temperature, velocity range and volume density, see Shetty et al. 2011a,b), but also on maximizing the dynamic range. This is critical for estimating the statistics of the underlying gas column density.

The PDFs of CO maps will be most informative to use when the optical depth of the tracer is known, as this will provide more certainty on their relationship with the parameters of turbulence. Conversely, if the sonic Mach number is known (and is supersonic), then one can get an order of magnitude estimate of the optical depth from the use of the statistics discussed in this work. Each tool has its own strengths and weaknesses. The Tsallis statistic is in general, more sensitive to the Alfvénic Mach number (especially in the case of strongly supersonic turbulence) than the moments and more sensitive to the changing radiative transfer parameter space. Ultimately, a use of multiple well-founded techniques, coupled with numerical simulations, will provide the most accurate parameters of the turbulence that exists in the observations. While PDFs were explored in this paper, they are not the only useful tool for studies of turbulence (see the Introduction). Therefore, in a subsequent paper we will extend the analysis to the impact of the radiative transfer on the spatial structure and spatial correlations.

Another question raised by this study is: could the PDF methods discussed in this work be used to find the observational sonic and Alfvénic Mach numbers, given that the optical depth is known? To answer this, we generate a new synthetic observation (with numerical resolution of ) as our test case. It has parameters =2.8, =2.0, and which are different from all the other models listed in Table 1 with which we used to draw our conclusions. However, if these parameters are now unknown, and all the information we have is the Figures of this paper, the integrated intensity map, and the fact that this case has an optical depth , what can we learn?

First, we can easily compute the statistical moments of the linear and logarithmic distributions. We find that the linear variance, skewness and kurtosis are 0.11, , and . The logarithmic variance is . Using the calculated values for the variance along with Equations 8 and Equations 9 with A=0.07 as given for the norm case (which has a slightly higher optical depth but is the closes in terms of the optical depth parameter space) in Table 2 provides a value of =3.0. Using the linear fit for kurtosis given in Burkhart et al. 2010 and shown in Figure 7 we find =2.7. As this is very close to the expected sonic Mach number, this again shows that the moments of cases with follow closely the statistics of the true column density. The skewness of the orange lines in Figure 6, also shows values consistent with the ranges given by skewness. For the estimation of the Alfvénic Mach number, we found that only incremental PDFs (i.e. the Tsallis fit) are substantially sensitive to , whereas the fully PDFs of the integrated intensity maps had primary sensitivities to the radiative transfer parameter space and sonic Mach number. We feel that multiple techniques are required in order to obtain both sonic and Alfvénic Mach numbers. This is especially true in the case of the Alfvénic Mach number, as many statistics have a primary sensitivity to the compressibility with the magnetic field contributing as a secondary influence. For example, future works should test how the velocity centroid anisotropy (see Esquivel & Lazarian 2011) is affected by radiative transfer.

This test case demonstrates that it is possible to determine parameters of turbulence with the use of statistics on CO data given information about the optical depth. In our case, we also had information about the driving of the turbulence in order to use parameter . Studies with different driving (i.e. b ranging from 1/3 to unity) should be considered in the future to assess the sensitivity of the present results on the driving mode of the turbulence. We stress that the reverse case of knowing the sonic Mach number, which can be provided by velocity dispersions (e.g. Kainulainen & Tan 2012), one can obtain information on the optical depth via use of PDFs. Ultimately, the more statistics that are used together, the more accurate the determination of the parameters of turbulence and radiative transfer environment will be. On this point, we note that PDFs represent only a small subset of statistics one can use to study turbulence with.

## 9 Conclusions and Summary

We investigate the PDF moments and Tsallis statistics of MHD turbulence including radiative transfer effects of the CO 2-1 transition in the context of their sensitivity to the sonic and Alfvénic Mach numbers and the radiative transfer environment (i.e. varying density and abundance values). We find that PDF descriptors applied to the integrated intensity maps generally show a large dependency on the radiative transfer parameter space in addition to the Mach numbers. In particular, we find that:

• PDFs of CO integrated intensity maps are generally log-normal, although deviations are present.

• The linear and logarithmic variance- relationship is observed in integrated intensity maps with optical depths an order of magnitude one or less. For cases with high CO abundances and densities and therefore optically thick lines, the variance looses sensitivity to the sonic Mach number. We recommend using the logarithmic variance for estimation of the sonic Mach number, as it is less sensitive to scaling (i.e. choice of X factor) than the linear variance.

• The higher order moments (skewness and kurtosis) of the integrated intensity distribution depend on the radiative transfer environment and sonic Mach number. For cases where the CO abundance and density values are too low to excite the full dynamic range of the gas, the skewness and kurtosis increase with increasing sonic Mach number more rapidly than the true column density. For very optically thick cases, the higher order moments lose their sensitivity to the sonic Mach number.

• The Tsallis fit parameters also show sensitivity to the radiative transfer environment however, dependencies on both sonic and Alfvenic Mach number exist in cases with an order of magnitude one or less. The Tsallis parameters and fit follow the true column density best for the case.

B.B. acknowledges support from the NSF Graduate Research Fellowship and the NASA Wisconsin Space Grant Institution. A.L. thanks both NSF AST 0808118 and the Center for Magnetic Self-Organization in Astrophysical and Laboratory Plasmas for financial support. All authors acknowledge support through grant SFB956/DFG and by the Deutsche Forschungsgemeinschaft, DFG, project number Os 177/2-1. This work was completed during the stay of A.L. as Alexander-von-Humboldt-Preisträger at the Ruhr-University Bochum and the University of Cologne.

## 10 Appendix

As mentioned in Section 3 the original simulations set the sound speed by varying the pressure and holding the RMS velocity constant. However, now we set a constant temperature of 10K for all the simulations. In this case, it makes sense to hold the sound speed constant and vary the RMS velocity in order to set the isothermal sonic Mach number. Here we rescale the velocity field of all the simulations in order to hold the sound speed constant for a given sonic Mach number. Changing the velocity field changes the overall radiative transfer calculation, however we find that the resulting PDFs and their dependencies on the sonic Mach number and optical depth are not affected. This is because the PDFs of the integrated intensity maps are not sensitive to velocity information, and are primarily sensitive to the optical depth and sonic Mach number. We demonstrate this in Figure 11, which is similar to Figure 3 except here we use simulations that vary the sonic Mach number via scaling the velocity field rather than the sound speed. We use the same radiative transfer parameter space (i.e. the same values of density and abundance) that was used in Figure 3 , however, because we alter the velocity field, the radiative transfer problem changes as well. Nevertheless, the behavior of the PDFs as a function of sonic Mach number and optical depth is unchanged compared with Figure 3. In both cases, the synthetic observations with optical depth on order unity display the most similarity to the column density PDF, while the cases with very low optical depth have PDF widths that are wider than the column density.

## References

• Armstrong et al. (1995) Armstrong, J. W., Rickett, B. J., & Spangler, S. R., 1995, ApJ, 443, 20
• Ballesteros-Paredes et al. (2007) Ballesteros-Paredes, J., Klessen, R. S., Mac Low, M.-M., & Vazquez-Semadeni, E., 2007, Protostars and Planets V, 63
• Bertoldi & Mckee (1992) Bertoldi & McKee, C., 1992, ApJ, 395, 140
• Burla & Vinas (2004a) Burlaga, L. F., & Vinas, A.-F., 2004a, Geophys. Res. Lett., 31, 16807
• Burla & Vinas (2004b) Burlaga, L. F., & F.-Vinas, A. 2004b, Journal of Geophysical Research (Space Physics), 109, 12107
• Burgh et al. (2007) Burgh, E., France, K., McCandliss S., 2007, ApJ, 658, 446
• Burla & Vinas (2005a) Burlaga, L. F., & -Vinas, A. F. 2005a, Physica A Statistical Mechanics and its Applications, 356, 375
• Burla & Vinas (2005b) Burlaga, L. F., & Vinas, A.-F. 2005b, Journal of Geophysical Research (Space Physics), 110, 7110
• Burla & Vinas (2006) Burlaga, L. F., & F.-Vinas, A. 2006, Physica A Statistical Mechanics and its Applications, 361, 173
• Berkhuijsen & Fletcher (2008) Berkhuijsen E., Fletcher, A., 2008, MNRAS, 390, 19
• Brunt (2010) Brunt, C. M., 2010, A&A 513, A67
• Brunt, Federrath, & Price (2010) Brunt, C. M., Federrath, C., & Price, D. J. 2010, MNRAS, 403, 1507
• Burkhart et al. (2009) Burkhart, B., Falceta-Goncalves, D., Kowal, G., Lazarian, A., 2009, ApJ, 693, 250
• Burkhart et al. (2010) Burkhart, B., Stanimirovic, S., Lazarian, A., Grzegorz, K., 2010, ApJ, 708, 1204
• Burkhart, Lazarian, & Gaensler (2012) Burkhart, B., Lazarian, A., Gaensler, B., 2012a, ApJ, 708, 1204
• Burkhart et al. (2012) Burkhart, B., Goodman, A., Lazarian, A., Rosolowski, E., 2012b, ApJ, submitted.
• Burkhart & Lazarian (2012) Burkhart, B., & Lazarian, A., 2012, ApJL, 755, 19
• Chepurnov et al. (2009) Chepurnov, A., Gordon, J., Lazarian, A., Stanimirovic, S., 2008, ApJ, 688, 1021
• Chepurnov & Lazarian (2009) Chepurnov, A., & Lazarian, A. 2009, ApJ, 693, 1074
• Chepurnov & Lazarian (2010) Chepurnov, A., & Lazarian, A., 2010, ApJ, 710, 853
• Cho, Lazarian, & Vishniac (2002) Cho, J., Lazarian, A., & Vishniac, E., 2002, ApJ, 566, 49
• Cho & Lazarian (2003) Cho, J. & Lazarian, A. 2003, MNRAS, 345, 325
• Elmegreen (1993) Elmegreen B. G., 1993, ApJ, 419, 29
• Elmegreen & Scalo (2004) Elmegreen, B. G. & Scalo, J.,2004, ARA&A,42,211
• Esquivel & Lazarian (2011) Esquivel, A., Lazarian, A., 2011, ApJ, 740, 117
• Eyink, Lazarian, & Vishniac (2011) Eyink, G., Lazarian, A., Vishniac E., 2011, ApJ, 743, 51
• Fed. et al. (2008) Federrath, C., Klessen, R. S., & Schmidt, W , 2008, ApJ, 688, 79
• Fed. et al. (2010) Federrath, et al., 2010, A&A, 512, A81
• Federrath & Klessen (2012) Federrath, C., & Klessen, R., 2012, ApJ, 761, 156
• Federrath & Klessen (2012) Federrath, C., & Klessen R., 2013, ApJ, 763, 51
• Gaensler et al. (2011) Gaensler et al. 2011, Nature, 478, 214-217
• Glover et al. (2010) Glover et al, 2010, MNRAS, 404, 2
• Goldreich & Sridhar (1995) Goldreich, P., & Sridhar, S., 1995, ApJ, 438, 763
• Goodman et al. (2009) Goodman, A. A., Rosolowsky, E. W., Borkin, M. A., Foster, J. B., Halle, M., Kauffmann, J., & Pineda, J. E. 2009a, Nature, 457, 63
• Goodman et al. (2009) Goodman, A. A., Pineda, J. E., & Schnee, S. L. 2009b, ApJ, 692, 91
• Hennebelle (2008) Hennebelle, P., & Chabrier, G. 2008, ApJ, 684, 395
• Hennebelle (2011) Hennebelle, P., & Chabrier, G., 2011, ApJ, 743, 29
• Heyer & Schloerb (1997) Heyer, M. H., & Schloerb, F. P. 1997, ApJ, 475, 173
• Heyer & Brunt (2004) Heyer, M., Brunt, C., 2004, ApJ, 615, 45
• Hill et al. (2008) Hill, A. et al., 2008, ApJ, 686, 363
• Hunter (1977) Hunter, C., 1977, ApJ, 218, 834
• Jeans (1902) Jeans, J. H., 1902, Phil. Trans. A., 199, 1
• Klessen et al. (2000) Klessen et al. 2000, ApJ, 535, 887
• Kainulainen & Tan (2012) Kainulainen, J., & Tan J., 2012, A&A, 549, A53
• Konstandin et al. (2012) Konstandin et al., 2012, ApJ, 761, 149
• Kowal et al. (2007) Kowal, G., Lazarian, A. & Beresnyak, A., 2007, ApJ, 658, 423
• Kowal & Lazarian (2010) Kowal, G., Lazarian, A., 2010, ApJ, 720, 742
• Krumholz & McKee (2005) Krumholz, M., McKee, C., 2005, ApJ, 630, 250
• Larson (1981) Larson, R. B., 1981, MNRAS, 194, 809
• Lazarian (2005) Lazarian, A., 2005, AIPC, 784, 42
• Lazarian (2006) Lazarian, A., 2006, ApJ, 645, 25
• Lazarian (2009) Lazarian, A., 2009, Space Sci.Rev., 143, 357
• Lazarian & Vishniac (1999) Lazarian, A., Vishniac, E., 1999, ApJ, 517, 700
• Lazarian, Santos-Lima,& de Gouveia Dal Pino (2010) Lazarian, A., Santos-Lima, R., & de Gouveia Dal Pino, E. 2010, Numerical Modeling of Space Plasma Flows, Astronum-2009, 429, 113
• Mckee & Ostriker (2007) McKee, C., & Ostriker E., 2007, ARAA, 45, 565
• Nordlund (1999) Nordlund, A., & Padoan, P. 1999, in Interstellar Turbulence, ed. J. Franco & A. Carraminana (Cambridge: Cambridge Univ. Press), 218
• Molina et al. (2012) Molina et al., 2012, MNRAS, 423, 2680
• Ossenkopf (2002) Ossenkopf, V., 2002, A&A, 391, 295
• Ossenkopf (2002) Ossenkopf, V., & Mac Low, M.-M., 2002, A&A, 390, 307 (Oss02)
• Ossenkopf, Krips, & Stutzki (2008) Ossenkopf, V., Krips, M., & Stutzki, J. 2008a, A&A, 485, 719
• Ossenkopf, Krips, & Stutzki (2008) Ossenkopf, V., Krips, M., & Stutzki, J. 2008b, A&A, 485, 917
• Ostriker, Stone, & Gammie (2001) Ostriker, E. C., Stone, J. M., & Gammie, C. F. 2001, ApJ, 546, 980
• Padoan & Jones (1997) Padoan, P. Nordlund, A., Jones, B. J. T., 1997, MNRAS, 288, 145
• Padoan & Nordlund (1999) Padoan, P., & Nordlund, A, 1999, ApJ,526, 279
• Padoan et al. (2007) Padoan, P., Nordlund, Å., Kritsuk, A. G., Norman, M. L., & Li, P. S. 2007, ApJ, 661, 972
• Price (2011) Price D., Federrath C., & Brunt C., 2011, ApJ, 727, 21
• Roman-Duval et al. (2011) Roman-Duval et al., 2011, ApJ, 740, 120
• Schneider et al. (2011) Schneider, N. et al., 2011, A&A, 592, A1
• Santos-Lima et al. (2010) Santos-Lima, R., Lazarian, A., de Gouveia Dal Pino, E. M., & Cho,J. 2010, ApJ, 714, 44
• Santos-Lima, de Gouveia Dal Pino, & Lazarian (2011) Santos-Lima, R., de Gouveia Dal Pino, E. M., & Lazarian, A. 2011,arXiv:1109.371
• Shetty et al. (2011a) Shetty, R.; Glover, S.C.; Dullemond, C.P.; Klessen R.; et al. 2011a, MNRAS 412, 1686
• Shetty et al. (2011b) Shetty, R.; Glover, S.C.; Dullemond, C.P.; Ostriker, E.C; et al. 2011b, MNRAS 415, 3253,
• Stutzki & Guesten (1990) Stutzki, J., & Guesten, R., 1990, ApJ, 356, 513
• Stutzki et al. (1998) Stutzki J., Bensch, F., Heithausen, A., Ossenkopf, V., & M. Zielinsky, M., 1998, A&A, .336, 697
• Spitzer (1979) Spitzer, L., Physical Processes in the Interstellar Medium, 1978, New York: John Wiley & Sons.
• Tofflemire, Burkhart, & Lazarian (2011) Tofflemire, B., Burkhart, B., Lazarian, A., 2011, ApJ,. 736,60
• Tsallis (1988) Tsallis, C., 1988,JSP, 52, 479
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters