{'fraction_non_alphanumeric': 0.033499611001784814, 'fraction_numerical': 0.03185208914923802, 'mean_word_length': 4.356416227478857, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Centroid Velocity Statistics of Molecular Clouds

We compute structure functions and Fourier spectra of 2D centroid velocity (CV) maps in order to study the gas dynamics of typical molecular clouds (MCs) in numerical simulations. We account for a simplified treatment of time-dependent chemistry and the non-isothermal nature of the gas and use a 3D radiative transfer tool to model the CO line emission in a post-processing step. We perform simulations using three different initial mean number densities of n_0 = 30, 100 and 300 cm^{-3} to span a range of typical values for dense gas clouds in the solar neighbourhood. We compute slopes of the centroid velocity increment structure functions (CVISF) and of Fourier spectra for different chemical components: the total density, H2 number density, 12CO number density as well as the integrated intensity of 12CO (J=1-0) and 13CO (J=1-0). We show that optical depth effects can significantly affect the slopes derived for the CVISF, which also leads to different scaling properties for the Fourier spectra. The slopes of CVISF and Fourier spectra for H2 are significantly steeper than those for the different CO tracers, independent of the density and the numerical resolution. This is due to the larger space-filling factor of H2 as it is better able to self-shield in diffuse regions, leading to a larger fractal co-dimension compared to CO.


INTRODUCTION
Understanding the dynamics of molecular clouds (MCs) is an important task in astrophysics. It is still a great debate which physical processes regulate star formation. Beyond self-gravity, the radiation field, magnetic fields and a complex thermodynamic and chemical evolution, several studies in the past have shown that turbulent motions play an important role for controlling the star formation process (see, e.g. Mac Low & Klessen 2004;Scalo & Elmegreen 2004;Elmegreen & Scalo 2004;McKee & Ostriker 2007;Hennebelle & Falgarone 2012).
Characterising interstellar turbulence is therefore important for understanding star formation. During the last years large improvements were made in understanding turbulence theory. Kolmogorov (1941) proposed a constant mean energy dissipation rate within an inertial range and assumed that all statistical averaged quantities only depend on this rate, leading to the law of Kolmogorov for incompressible turbulence. This directly leads to the famous 5/3-law of the energy spectrum with a statistically quasi-equilibrium state. Nevertheless, many experimental and theoretical studies in the past decades have shown that the assumption of a constant mean energy dissipation rate is not accurate (Benzi et al. 1993;She & Leveque 1994;Dubrulle 1994;Frisch & Kolmogorov 1995) and seems to deviate significantly from measurements for higher orders of the structure function. This discrepancy has been referred to as intermittency corrections to the Kolmogorov theory. A number of phenomenological models have been proposed since then to explain the differences. She & Leveque (1994) proposed an intermittent scaling for filamentary structures, in excellent agreement with tunnel-flow experimental data (Benzi et al. 1993) for incompressible turbulence. Boldyrev et al. (2002) expanded the work of She & Leveque (1994) for sheet-like shocks or a general fractal dimension of the dissipative structures for compressible turbulence. Later on, other studies developed models for a better understanding of turbulence involving magnetic fields (see, e.g. Goldreich & Sridhar 1995;Cho et al. 2002).
Studying interstellar turbulence is difficult, since observations are always a complex convolution of the density and the velocity field, affected by many other astrophysical processes. Assessing the influence of projection effects relies on numerical simulations of turbulence that can be compared to observational measurements. In that sense, Lis et al. (1996) have shown that it is possible to use centroid velocity increments (CVI) in order to trace the most intense velocity structures of regions in the interstellar medium (ISM), which do not form stars and to analyse intermittency effects on the basis of the two-point statistics.
Recent studies have revealed in great detail that centroid velocities (CV) and their increments are sensitive to the turbulent driving, density fluctuations and temperatures along the line-of-sight (Lazarian & Esquivel 2003;Ossenkopf et al. 2006;Esquivel et al. 2007;Hily-Blant et al. 2008;Federrath et al. 2010). Furthermore, the non-uniform chemical makeup of the gas and optical depth effects might also alter the CVI statistics (see, e.g. Burkhart et al. 2013b andBurkhart et al. 2014). Chemical effects may be important for an accurate determination of the velocity fluctuations, since the relation between H2 as the main constituent of the ISM and CO as a tracer is not linear (see, e.g. Glover et al. 2010;Shetty et al. 2011a). Optical depth effects may as well affect projected structures due to different opacities of the gas tracers (see, e.g. Burkhart et al. 2013a,b, 2014and Lazarian & Pogosyan 2004. Nevertheless, the question of how different physical processes affect the overall two-point statistics of the CVI remains elusive although much progress has been made during the last years. In this paper we analyse high-resolution, threedimensional chemistry models in numerical simulations of turbulent interstellar gas (Glover et al. 2010) to get a better insight into how the choice of chemical species as a tracer of gas affects the two-point statistics of the CV. Therefore we use time-dependent chemistry for modelling the formation and destruction of molecules in the ISM. We then perform radiative transfer calculations and convert our simulations into synthetic 2D CV maps. We focus in this paper on how chemistry and radiative transfer affect the CVI turbulence statistics and the Fourier spectra in different environments, while spanning a range of possible densities representative of MCs.
In Sec. 2 we present the simulations, the methods used to compute the CVI, the Fourier spectra and the radiative transfer post-processing. In Sec. 3 we present the results of our studies with different mean number densities. We discuss our results in Sec. 4 and present our conclusions in Sec. 5.

MHD and chemistry simulations
The simulations that we examine in this paper were performed using a modified version of the Zeus-Mp magnetohydrodynamical code (Norman 2000;Hayes et al. 2006). We make use of a detailed atomic and molecular cooling function, described in detail in Glover et al. (2010) and Glover & Clark (2012), and a simplified treatment of the molecular chemistry of the gas. Our chemical treatment is based on the work of Nelson & Langer (1999) and Glover & Mac Low (2007), and allows us to follow the formation and destruction of H2 and CO self-consistently within our simulations. Full details of the chemical model can be found in Glover & Clark (2012).
We consider a periodic volume, with side length 20 pc, filled with uniform density gas. We run simulations with three different initial densities, n0 = 30, 100 and 300 cm −3 , where n0 is the initial number density of hydrogen nuclei. The initial temperature of the gas is set to 60 K. We initialise the gas with a turbulent velocity field with uniform power between wavenumbers k = 1 and k = 2 (Mac Low et al. 1998;Mac Low 1999). The initial 3D rms velocity is vrms = 5 km s −1 , and the turbulence is driven so as to maintain vrms at approximately the same value throughout the run. The gas is magnetised, with an initially uniform magnetic field strength B0 = 5.85 µG, initially oriented parallel to the z-axis of the simulation. We neglect the effects of self-gravity. We assume that the gas has a uniform solar metallicity and adopt the same ratios of helium, carbon and oxygen to hydrogen as in Glover et al. (2010) and Glover & Mac Low (2011). At the start of the simulations, the hydrogen, helium and oxygen are in atomic form, while the carbon is assumed to be in singly ionised form, as C + . We also adopt the standard local value for the dust-to-gas ratio of 1:100 (Glover et al. 2010), and assume that the dust properties do not vary with the gas density. The cosmic ray ionisation rate was set to ζ = 10 −17 s −1 , but we do not expect our results to be sensitive to this choice (Glover et al. 2010). For the incident ultraviolet radiation field, we adopt the standard parameterisation of Draine (1978). More detailed information about the simulations that we examine in this paper can be found in Bertram et al. (2014). The runs analysed in this paper were all performed with a numerical resolution of 512 3 zones. We discuss the influence of our limited numerical resolution in more detail in Appendix A.

Radiative transfer
To model the CO (J = 1 → 0) line in the MCs for both 12 CO and 13 CO we employ the radiative transfer code Radmc-3d 1 (Dullemond 2012). We use the Large Velocity Gradient (LVG) approximation (Sobolev 1957) to compute the level populations. The LVG implementation in Radmc-3d is described in Shetty et al. (2011a). For more information about the usage of Radmc-3d on our data we refer the reader to Bertram et al. (2014).
Our simulations do not explicitly track the abundance of 13 CO, and so we need a procedure to relate the 13 CO number density to that of 12 CO. A common assumption is that the ratio of 12 CO to 13 CO is identical to the elemental abundance ratio of 12 C to 13 C (see, e.g. Roman-Duval et al. 2010). In the majority of our analysis, we make the same assumption and set the 12 CO to 13 CO ratio to a constant value, R 12/13 = 50.
However, strictly speaking, we expect the 12 CO to 13 CO ratio to faithfully reflect the elemental 12 C to 13 C ratio only in gas where essentially all of the carbon is incorporated into CO. When this is not the case (e.g. at low density and/or low extinction), the effects of chemical fractionation (Watson et al. 1976) and selective photodissociation of 13 CO (see, e.g. Visser et al. 2009) can significantly alter the 12 CO/ 13 CO ratio (Röllig & Ossenkopf 2013;Szűcs et al. 2014). To explore the effect that this may have on the statistical properties of the 13 CO line emission, we have produced a set of 13 CO emission maps, where we attempt to account for spatial variations in the 12 CO/ 13 CO ratio. To do this, we made use of the recent numerical results of Szűcs et al. (2014). They find that in turbulent molecular clouds, there is a tight relationship between the mean value of R 12/13 along a given line-of-sight and the 12 CO column density along the same line-of-sight, and give a simple polynomial fitting formula for R 12/13 as a function of the 12 CO column density. To produce our second set of 13 CO maps, we therefore first calculate the 12 CO column density along one axis of our datacubes, and then use these values to compute the mean value of R 12/13 for each line-of-sight using the Szűcs et al. (2014) fitting formula. Finally, we derive the number density of 13 CO in each cell along each line-of-sight by dividing the 12 CO number density by the mean value of R 12/13 . These 13 CO number densities are then used to compute the 13 CO emission in the same fashion as in our constant ratio models. We note that although this procedure neglects any spatial variations in R 12/13 along a given line-of-sight, Szűcs et al. (2014) have shown that the resulting emission maps differ by only a few percent from those derived from a self-consistent (but costly) chemical treatment. The differences in the CVI statistics of the two sets of 13 CO emission maps -those derived using a constant R 12/13 and those that use a spatially varying value of R 12/13 -will be discussed in Sec. 3.3 below.
Although Radmc-3d accounts for the dust continuum emission as well as the CO line emission, we make sure to subtract off the continuum term before analysing the statistical properties of the synthetic maps. For further information about the usage of Radmc-3d and the calculation of the dust continuum and the line emission we refer the reader to Dullemond (2012).
The radiative transfer calculation yields positionposition-velocity (PPV) cubes of brightness temperatures TB, which are used to compute centroid velocities. We will refer to these as the "intensity" models. In analogy, we construct centroid velocity maps out of the PPP simulation data of the density and velocity field, as described in Sec. 2.3. We will refer to these as the "density" models.

Centroid velocity increments
Our study is based on the two-point statistics of line centroid velocities. Following Lis et al. (1996), we analyse the centroid C(x, y) = C(r) of line-of-sight projected velocities, defined as where the variable vz(r, z) is the line-of-sight velocity in the z-direction. The quantity F (r, z) is a statistical weight and denotes either the underlying density field or the brightness temperatures from the radiative transfer post-processing. The centroid velocity increment δC is defined as the separation between two points of the CV by spatial distances , i.e.
It thus connects centroid velocities of different regions in the plane.
There are many other techniques to measure velocity fluctuations as a function of spatial scale, e.g. the spectral correlation function (Rosolowsky et al. 1999), the velocity channel analysis (Lazarian & Pogosyan 2000, 2004, the ∆variance method (Ossenkopf et al. 2008a,b)   . CVI structure function for orders up to p = 6 with the inferred slopes ζp as a function of spatial scale for one snapshot in time, normalised by the total box size. The velocity here is weighted by the total density field for an initial number density of n 0 = 300 cm −3 . The fitting range is indicated by a horizontal dashed line. Error bars denote 1-σ variations for all 3 line-of-sight directions.
Heyer 2002). The latter has been applied to study the gas dynamics of the same simulation data used in this paper in Bertram et al. (2014).

Structure functions of the CVI
Statistical moments p of the distribution of all possible CVI are called structure functions (CVISF). They are defined as where the average is computed over all possible CVI separated by in the plane. When p increases, the CVISF give a stronger weight to rare events. In our study we have computed the CVISF for each possible line-of-sight direction x, y and z. Although we use a weak magnetic field in the zdirection, the turbulence in our simulations is trans-Alfvénic or mildly super-Alfvénic (see Table 1). The field lines are therefore dragged along with the turbulent flow, with the result that the turbulence remains approximately isotropic (see, e.g. Burkhart et al. 2014). Hence, we do not expect significant variations of the CVISF along the different linesof-sight. We therefore average all CVISF of the different lineof-sight directions and compute all CVISF for p 6. Finally, the CVISF were fit to power laws of the form where ζp is a function that depends on the statistical moment p of the CVI distribution. To calculate the scaling exponents, we use a fitting range from 1/16 to 1/5 of the total box size, as used by Federrath et al. (2010) and . For a box with D = 512 grid cells for each side this translates to 32 and 102 cells. Figure 1 shows an example of a CVISF for orders up to p = 6 for one snapshot in time, where the velocity is weighted by the total density field from model n300.

Fourier and power spectra
Fourier spectra are used to analyse the scale dependence of the kinetic energy of turbulence simulations. Assessing 3D turbulence power spectra, Kolmogorov (1941) showed through statistical arguments that the kinetic energy should scale as E(k) ∝ k −2/3 in case of incompressible turbulence, where k is the wavenumber. Later, Burgers (1948) derived E(k) ∝ k −1 for supersonic, compressive turbulence. We define the power spectrum as the change of the kinetic energy with its wavenumber, i.e. P (k) = dE/dk. The kinetic energy spectrum can be computed in Fourier space via where v denotes the Fourier transform of the velocity field and v * its complex conjugate. Following Kolmogorov (1941), the kinetic energy down the cascade in the inertial range scales with δv 2 ∝ ζ 2 , where ζ2 is the slope of the second order structure function obtained in the inertial range (or fitting range, as we term it in our case). Thus, the power spectrum can be computed and we obtain where we have used k ∝ −1 and set α = −1−ζ2 as the value of the spectral slope. Hence, the power spectrum can be determined by evaluating ζ2 of the second order structure function. Using Eq. 6, we then obtain slopes P (k) ∝ k −5/3 and P (k) ∝ k −2 for both Kolmogorov (α = −5/3) and Burgers (α = −2) turbulence, respectively. In this paper, we compute power spectra for 2D maps of CV in order to obtain information about the underlying physics of the different chemical components.

Theoretical relations for the SF slopes
Following Kolmogorov (1941), we assume that the kinetic energy scales as E(k) ∝ δv 2 and E(k) ∝ k −2/3 in the case of incompressible turbulence, from which it directly follows that δv ∝ 1/3 . Furthermore, by taking δv p ∝ ζp , we can infer the theoretical scaling relation for the structure functions (SF) in the sense of Kolmogorov and find This law will be referred to as K41. Later on, She & Leveque (1994) applied intermittency corrections to K41 and set ζp = p/3 + τ p/3 , since Eq. 7 seems to deviate significantly from many experimental studies for p > 3 (Benzi et al. 1993;She & Leveque 1994;Dubrulle 1994;Frisch & Kolmogorov 1995). The function τ p/3 thereby includes additional terms that depend on p. She & Leveque (1994) propose which we refer to as the SL94 model. As shown by Boldyrev et al. (2002), hereafter B02, and Dubrulle (1994), one important parameter in models of supersonic turbulence is the fractal co-dimension, with D = 2. Hence, Boldyrev et al. (2002) and Dubrulle (1994) generalised Eq. 8 and found where Σ = 1−Θ/C. The other parameters Θ and ∆ enter the theory through the so-called Kolmogorov refined similarity hypothesis, relating the scaling of fluctuations of the velocity field to the scaling of the fluctuating energy dissipation. We leave those two as free arbitrary parameters for now.

RESULTS
We analyse different numerical models with and without radiative transfer post-processing and vary the initial number density in the simulation box, in order to study the influence on the structure functions. Table 1 gives an overview of the numerical models that we study in this paper. We quote both mass-and volume-weighted quantities, which we define using the expressions f mass = f ρdV / ρdV and f vol = f dV / dV , respectively. More detailed analyses of the H2 and CO distributions produced in this kind of turbulent simulation can be found in previous studies by Glover et al. (2010) and Shetty et al. (2011a). Examples of the integrated intensity and column density PDFs are shown in Shetty et al. (2011a). Figure 2 shows the centroid velocity maps computed via Eq. 1 for all models and for all chemical components.

Analysis of the CVISF
We compute the CVISF for all models (n30, n100, n300) and all chemical tracers: the total density, H2 density, CO Figure 2. Centroid velocities in z-direction for models n30 (left column), n100 (middle column) and n300 (right column). From top to bottom: CV maps of the total density, H 2 and 12 CO density followed by CV maps of the 12 CO and 13 CO integrated intensities. Each side has a length of 20 pc. Note that the velocity field of the n30 model uses a different turbulent seed than the n100 and the n300 model. Black areas in the 12 CO and 13 CO map of the n30 model denote regions where the brightness temperatures are zero along the LoS. density, 12 CO and 13 CO intensity. We average slopes of the CVISF using 3 snapshots in time (with 3 line-of-sight directions each) where we can assume both the chemistry and the turbulence to be in a stationary and converged state. In Figure 3 we show slopes ζp of the CVISF as a function of the order p for all models and chemical tracers. The error bars denote temporal and spatial 1-σ fluctuations from the different time snapshots. Moreover, we also show the K41, SL94 and B02 scaling relations in the plots. We have to keep in mind that a direct comparison of the theoretical models with our CV statistics is nontrivial and sometimes impossible, since we analyse velocity centroids, while the theory predicts scaling relations for the full turbulent velocity field. For a more detailed discussion about the influence of the projection, we refer the reader to Sec. 4.3 and 4.4.
We find differences between the various models and tracers. In our n300 model we obtain similar slopes of the total density and the H2 density model as well as for the CO density and the 13 CO intensity model within their error bars. The slopes of the total density and the H2 density model are steeper compared to the different CO tracer models, which flatten to higher order p of the CVISF. The slopes ζp of the 12 CO intensity model make up an intermediate case. They are slightly steeper than the CO density and 13 CO intensity model, but flatter than those of the total density and H2 density model. Regarding the n100 model, we find consistent values for ζp for all CO tracer cases and for all orders p. Furthermore, we again measure steeper slopes for the total density and H2 density model, which now slightly deviate from each other at higher order of the CVISF.
In our n30 model we obtain a larger discrepancy between the total density and the H2 density model. Although the temporal and spatial 1-σ fluctuations are large for both cases, we find a clear separation of the slopes for higher p. Again, all different CO tracers agree with each other within their error bars and are significantly flatter and those of the H2 and total density cases.

Analysis of the Fourier spectra
Following Eq. 5, we compute Fourier spectra and the scaling exponents of the power spectra by fitting a linear function in log-log space to the Fourier spectra in the fitting range of our simulations. The fitting range in k-space translates to wavenumbers from 5 to 16 and has to be chosen carefully due to the artificial accumulation of energy further down the cascade, known as the bottleneck effect (Falkovich 1994;Dobler et al. 2003;Kaneda et al. 2003;Haugen & Brandenburg 2004;Kritsuk et al. 2007;Beresnyak & Lazarian 2009). In analogy to Figure 3, Figure 4 shows with k 2 compensated Fourier spectra for the different density models and all chemical components.
In the n300 model we again find a good correlation between the total density and the H2 density. We measure slopes α of −1.62 ± 0.09 for the former and −1.60 ± 0.09 for the latter, which agree well within their errors. For the CO density, 12 CO and 13 CO intensity cases we find slopes of −1.27 ± 0.10, −1.36 ± 0.08 and −1.30 ± 0.09, which do not show significant differences within their errors.
The same trends are seen in the other density models. In the n100 model we again do not find a significant difference Figure 3. Slopes ζp of the structure functions against order p for all chemical models, i.e. the total density, H 2 density and the CO density model as well as the 12 CO and 13 CO intensity model. For comparison, we also show the theoretical scaling relations (black lines) from Kolmogorov (1941), She & Leveque (1994) and Boldyrev et al. (2002), denoted as K41, SL94 and B02. From top to bottom: runs of the different initial number densities, i.e. n 0 = 300, 100 and 30 cm −3 . Error bars denote temporal and spatial 1-σ fluctuations. Please note that a direct comparison of the theoretical models with our CV statistics is complicated due to the projection along the LoS. For a more detailed discussion about the influence of the projection, we refer the reader to Sec. 4.3 and 4.4. between the total density and the H2 density, for which we measure slopes α of −1.57 ± 0.09 and −1.54 ± 0.09, respectively. For the CO density, 12 CO and 13 CO intensity cases we find very flat slopes of −0.98 ± 0.07, −1.01 ± 0.08 and −0.87 ± 0.09.
In the low-density n30 model, the slopes α of the total density (−1.49 ± 0.11) and H2 density (−1.41 ± 0.11) cases agree well within their errors, while the slopes of the CO density, 12 CO and 13 CO intensity cases are again much flatter and we measure −1.04±0.09, −1.01±0.10 and −1.00±0.10.
Overall, we find a significant difference in the slopes of the H2 density and all CO tracers for all three density models. Slopes of the total density and H2 density are found to be significantly steeper than those produced by the different CO tracers, while we do not obtain significant differences in the slopes between the different CO tracers within their errors.

Variation of the abundance of 13 CO
Our numerical simulations follow the chemistry of 12 CO, but not of its isotope 13 CO. Consequently, in order to generate the data on the number density of 13 CO that we need for our radiative transfer calculation, it is necessary to specify a conversion factor between 12 CO and 13 CO. As we explain in Sec. 2.2, in most of our analysis, we assume a fixed 12 CO to 13 CO ratio, R 12/13 = 50. However, in real molecular clouds, we expect spatial variations in R 12/13 in any regions in which the carbon is not all locked up in CO, caused by the two effects of chemical fractionation and selective photodissociation (Röllig & Ossenkopf 2013;Szűcs et al. 2014).
To explore the effect that these variations may have on the CVI statistics of the 13 CO emission maps, we produce a set of 13 CO emission maps, using a spatially varying R 12/13 generated following the prescription given in Sec. 2.2 and in Szűcs et al. (2014). After applying the radiative transfer post-processing to both models, we compute the two maps of the CV and compare the slopes of the CVISF.
Although we find slight differences in the CV maps on small scales, the slopes ζp of the CVISF are essentially the same in both cases. The differences between the results in the two models are smaller than the 1-σ temporal and spatial errors in the different snapshots. We thus conclude that the use of a constant 12 CO/ 13 CO ratio is sufficient for obtaining proper slopes of the CVISF for 13 CO. Figure 3, we find large differences in the CV-ISF slopes for the different chemical components and density models. In the n300 model, we find a good correlation between the total density and the H2 density within their error bars. This is because about 98% of the initial atomic hydrogen gas is in molecular form at this time, meaning that H2 is an excellent tracer of the total density. Furthermore, 13 CO is optically thin and thus we expect the 13 CO intensity PPV model to trace the underlying 3D PPP CO density model well, which we can see in Figure 3. This argument is also valid for the other density models. On the other hand, Figure 4. Fourier energy spectra multiplied by k 2 as a function of scale k for all density models and chemical components, i.e. the total density, H 2 density and the CO density model as well as the 12 CO and 13 CO intensities. From top to bottom: runs of the different initial number densities, i.e. n 0 = 300, 100 and 30 cm −3 . Error bars denote temporal and spatial 1-σ fluctuations. In all models, the total density and H 2 density cases show a significantly steeper slope compared to the CO tracer components. The fitting range is indicated by a horizontal dashed line. Slopes with the fitting errors are given in each plot for the different species.
12 CO is optically thick and traces H2 well only in dense regions (see Bertram et al. 2014). This is why it lies between the total density and the CO density models. We discuss the physical implications for this result in Sec. 4.4. Turning to the low-density model n30, we find a larger discrepancy in the slopes of the H2 density and the total density field, because only 61% of the initial atomic hydrogen gas is now in molecular form. Thus, the H2 gas does not trace the total gas very well, although it is still the best tracer for the total density field in general. Looking at the different CO tracers, we find a good correlation between the CO density, 12 CO and 13 CO intensity cases within the error bars.
The n100 model is intermediate between n30 and n300. We find a better correlation of the H2 and the total density cases than in the n30 model, which comes from the higher fractional abundance of H2 compared to the n30 model. Furthermore, we again obtain a good correlation between all CO tracer components within their error bars, which means that the turbulence statistics in all three cases are similar.
Nevertheless, regarding the question of whether CO tracer molecules can be used to infer the dynamics of the underlying total density or H2 density field, we find that in all density models, the slopes of the total density or the H2 density field might be underestimated by a factor up to 3 by the use of different CO tracers.

Slopes of the Fourier spectra
We can understand the significant differences between the H2 and the CO tracer cases in the Fourier spectra of the distinct density models by comparing Figure 3 with Figure 4. In Eq. 6 we have shown that the slopes α of the Fourier spectra are related to the structure functions of the turbulence over the second order structure function ζ2 via α = −1 − ζ2. The interpretation of slope differences ∆α = ζ in Figure  4 between different chemical tracers is thus related to the differences in slopes of the second order structure functions. If the differences of values for ζ2 within one density model are large, we also expect large differences to occur in the slopes of the Fourier spectra, which is indeed the case. The "average" difference for ζ2 between H2 and the CO tracers is largest in the n100 model. Hence, we also find the largest slope differences between those two cases in our Fourier spectra for the n100 model. This is why Fourier slopes in Figure  4 of all CO tracers are significantly flatter then those of H2. Since we always find the highest slope values for the total density and H2 density in Figure 3, we also measure steeper exponents of those components in the Fourier spectra.
However, we have to keep in mind that we have used 3D statistical measures to understand the behaviour of the 2D centroid velocity maps. Nevertheless, since we only want to compare relative scaling behaviours in our models, we can safely use 3D statistical tools of turbulence analysis in order to work out trends in 2D data available for observers.
Furthermore, Lazarian & Pogosyan (2004) predicted that the spectral slope should saturate to a value of -3 for optically thick data in the context of integrated intensities, which was confirmed by several other studies in the past (see, e.g. Stutzki et al. 1998;Stanimirović & Lazarian 2001;Lazarian & Pogosyan 2006;Padoan et al. 2006;Burkhart et al. 2010, Burkhart et al. 2013a. We note that in the context of Lazarian & Pogosyan (2004), the power spectra are defined via E(k) ∝ k 2 P (k) ∝ k β . Thus, if we want to compare our results to those found in Lazarian & Pogosyan (2004), we have to subtract a value of -2 in order to compensate for the k 2 scaling, i.e. β = α − 2, where α denotes the Fourier slopes given in Figure 4 and β the slopes listed in Lazarian & Pogosyan (2004) or Burkhart et al. 2013a. Table  2 gives an overview about the spectral slopes β computed from our slope values α given in Figure 4. Independent of the individual component and the initial density, we also observe a saturation of the spectral slope indices to a value of β ≈ −3 for CO (see also Lazarian & Pogosyan 2004). However, we have to keep in mind that our statistics deal with centroid velocities, while previous studies elaborated the dependence of the spectral slopes on integrated intensities. Thus, we expect some modifications as in our case there is a dependency on both density and velocity, which could explain the slight steepening of slope values to higher densities seen in Table 2.

Projection effects
A direct comparison of CVISF to models of 3D turbulence such as proposed by She & Leveque (1994) or Boldyrev et al. (2002) is complicated, since maps of centroid velocities are a complex convolution of the density (or brightness temperature) with the velocity field. Moreover, previous studies have shown that physical effects like opacity or the sonic Mach number also significantly influence the result of the CV analysis (see, e.g. Lazarian & Pogosyan 2004, Burkhart et al. 2013aand Burkhart et al. 2014). Nevertheless, CVI statistics provide a useful tool to study effects of turbulence in observational measurements and to compare those to numerical experiments, since observations rely on spectral PPV information from which the CV maps can directly be inferred.
As shown by Lazarian & Esquivel (2003), Ossenkopf et al. (2006), or Esquivel et al. (2007, CVI statistics are significantly influenced by the underlying density field and thus differ from the pure 3D velocity statistics, unless the ratio of the density dispersion to the mean density is small. This is usually the case in subsonic flows where the Mach number is supposed to be small, which is clearly not the case in our simulations, where we measure supersonic Mach numbers. Table 1 gives both the mean Mach number and the ratio of density dispersion and mean density, σρ/ ρ . We find values σρ/ ρ 3 and therefore expect the mutual convolution of density and velocity field to play a significant role. Ossenkopf et al. (2006) give an upper limit of σρ/ ρ 0.5 in order to get a proper matching of the CVISF with the 3D structure functions of the turbulence. We thus expect a significant loss in the signatures of intermittency in our CVI statistics compared to the proper 3D data, which limits a direct comparison of the CVISF to the different moments of the structure functions of the underlying 3D velocity field. For further information about the density variance-Mach number relation we refer the reader to Vázquez-Semadeni & García (2001), Lemaster & Stone (2008), Federrath et al. (2008, Molina et al. (2012), Burkhart & Lazarian (2012) and .

Connection to the fractal dimension
As shown in Figure 3, we find increasing differences between the H2 and the CO tracers with higher order p of the CVISF for all density models. Although we analyse 2D CVI statistics, we can nevertheless gain useful information about the most dissipative structures in our simulations using theoretical models of 3D velocity statistics. However, as discussed in Kowal et al. (2007), we have to keep in mind that supersonic turbulence contaminates the scaling relations of the CV statistics due to a complex convolution of the supersonic density with the velocity.
In the context of She & Leveque (1994) (C = 2), Eq. 9 predicts larger slope values ζp than in the model given by Boldyrev et al. (2002) (C = 1), in which ζp becomes flatter for higher order p. Thus, the increasing differences between the H2 and the CO tracers with higher order p of the CV-ISF seem to be related to the fractal co-dimension C, which can be interpreted as a space-filling factor, describing the distribution of gas in the MC. As argued in Bertram et al. (2014), H2 is more extended than CO as it is better able to self-shield in diffuse regions, while CO is a good tracer of H2 only in compact regions. Hence, we would predict a larger fractal co-dimension for the H2 gas than for the various CO tracers, leading to the increasing differences in ζp as observed in Figure 3 between H2 and CO. However, due to projection and resolution effects it is rather complicated to give general values for the fractal co-dimension of the 3D velocity statistics. The question of how the fractal dimensions of 2D and 3D velocity statistics are related with each other can be addressed in a follow-up study.

Limitations of the models
As described in Bertram et al. (2014), there are some assumptions and limitations inherent to our numerical models that we must keep in mind when interpreting the results. For example, our models do not account for self-gravity, star formation, stellar feedback (e.g. by SN, stellar radiation, etc.) or large-scale dynamics (e.g. spiral arms, SN shells, etc.). Furthermore, we find a minor resolution dependence of the CVISF slopes for CO as discussed in Appendix A, which limits the ability to directly compare our theoretical values with observational measurements. However, we can gain useful information about the underlying physics through the relative scaling of the CVISF slopes to get an idea of how the chemistry and a radiative transfer affect the trends seen in tools like CVI statistics or Fourier spectra, which are commonly used by observers.
For future investigations, higher resolution models would be needed in order to measure converged slope values that can be compared to observations. Additionally, simulations that span a wider range of physical parameters, i.e. different metallicities or radiation fields, are important to study in order to analyse the influence on the CVI statistics. Star formation, feedback and stellar winds could also help to improve our understanding of turbulence in MCs.

SUMMARY AND CONCLUSIONS
We analysed centroid velocity increment structure functions (CVISF) and Fourier spectra of MCs with time-dependent chemistry and a radiative transfer post-processing for models of different initial number densities and chemical components: the total number density, H2 and CO density (each without radiative transfer) as well as 12 CO and 12 CO intensity (both with radiative transfer). In each case, we computed CVISF and their Fourier spectra and analysed the slopes of the CVISF within a fitting range used by Federrath et al. (2010) and  in order to study the influence of the chemistry and the radiative transfer on the results of the CVI statistics. We report the following findings: • Optical depth effects can be important for the CVI structure analysis. We find different behaviours of the CV-ISF slopes for 12 CO and 13 CO in different density environments (see Sec. 4.1).
• We find the slopes of H2 in Fourier spectra generally to be steeper than the slopes of different CO tracers (see Sec. 4.2). We also find the slopes of the CVISF for the total density and H2 density to be steeper by a factor of up to 3 than the slopes of different CO tracers.
• We find β −3 for all spectral slope values measured in this study (see Sec. 4.2), saturating at a value of -3 in the synthetic CO emission data cubes (see also Lazarian & Pogosyan 2004).
• We expect the CO gas to have a significantly smaller fractal co-dimension than the H2 gas, which means that it is less space-filling than H2 (see Sec. 4.4).
• We do not find any variations in the slopes of the CVISF using more realistic abundances of 13 CO in MCs compared to a constant 12 CO/ 13 CO scaling ratio throughout the whole cloud (see Sec. 3.3).
• Following the results of previous studies (see Sec. 4.3), we expect a significant loss of information in the turbulence statistics due to projection, since we find values of σρ/ ρ 3 in our simulations of supersonic turbulence. This limits a direct test of the CVISF to the different moments of the 3D structure functions. Milky Way System", and the SPP (priority program) 1573, "Physics of the ISM". Most of the simulations analysed in this paper were performed using the Ranger cluster at the Texas Advanced Computing Center, using time allocated as part of Teragrid project TG-MCA99S024. Additional simulations were performed on the kolob cluster at the University of Heidelberg, which is funded in part by the DFG via Emmy-Noether grant BA 3706, and via a Frontier grant of Heidelberg University, sponsored by the German Excellence Initiative as well as the Baden-Württemberg Foundation. RSK acknowledges support from the European Re- To study the influence of resolution on the CVISF slopes and the Fourier spectra, we have performed runs of 512 3 and 256 3 grid cells and evaluated the CVISF and Fourier slopes for all available chemical components and initial number densities. As an example, Figure A1 shows slopes of the CVISF for an initial number density of n0 = 100 cm −3 for both a resolution of 512 3 and 256 3 grid cells. Table A1 gives the corresponding slope values and errors for both resolution models for p = 6, since we would expect variations of the CVISF slopes due to intermittency particularly at higher orders. For the total density and the H2 density, we measure similar slopes within the errors, while the slopes of the different CO tracers differ by 1-2σ. This sensitivity to numerical resolution is a consequence of the high degree of chemical inhomogeneity in the numerical simulations. As shown in Bertram et al. (2014), CO is mainly located in dense gas regions which can be resolved more accurate in the 512 3 than in the 256 3 runs. It is thus rather complicated to compare non-converged CVISF slopes with observational measurements. Nevertheless, we find that although the individual slope values might be different for the distinct resolutions, the relative scaling behaviour of ζp is conserved throughout the different resolution runs and hence our fundamental physical conclusions derived in this paper should be unaffected. The same arguments hold for the Fourier spectra given in Figure A2. Although the slope values of the individual chemical components are significantly different between the two resolution models, the relative scaling behaviour of the energy spectra is conserved. Figure A1. Slopes ζp of the structure functions against order p for all chemical models, i.e. the total density, H 2 density and the CO density model as well as the 12 CO and 13 CO intensity model. From top to bottom: runs of different simulations of 512 3 and 256 3 grid cells with a fixed initial number density of n 0 = 100 cm −3 . Error bars denote temporal and spatial 1-σ fluctuations. Although the individual slope values might be different for the distinct resolutions, the relative scaling behaviour of ζp is conserved.

Resolution
Total density H 2 density 12 CO density 12 CO intensity 13 CO intensity 512 3 1.66 ± 0.13 1.45 ± 0.15 0.83 ± 0.10 0.84 ± 0.13 0.88 ± 0.05 256 3 1.68 ± 0.10 1.45 ± 0.08 0.72 ± 0.06 0.75 ± 0.09 0.74 ± 0.07 Table A1. Slopes of the CVISF for our highest order p = 6 for the different chemical components and for different runs with 512 3 and 256 3 grid cells. As an example, we show values for a fixed initial number density of n 0 = 100 cm −3 . The slopes of the total density and H 2 density cases are similar within their errors, while the slope values of the different CO tracers significantly deviate with resolution. Figure A2. Fourier energy spectra multiplied by k 2 as a function of scale k for all chemical components, i.e. the total density, H 2 density and the CO density model as well as the 12 CO and 13 CO intensities. From top to bottom: runs of different simulations of 512 3 and 256 3 grid cells with a fixed initial number density of n 0 = 100 cm −3 . Error bars denote temporal and spatial 1-σ fluctuations. In all models, the total density and H 2 density cases show a significantly steeper slope compared to the CO tracer components. The fitting range is indicated by a horizontal dashed line. Slopes with the fitting errors are given in each plot for the different species. Although the individual slope values might be different for the distinct resolutions, the relative scaling behaviour of the energy spectra is conserved.