# Structure analysis of simulated molecular clouds with the $\Delta$-variance

Erik Bertram1, Ralf S. Klessen1,2,3 & Simon C. O. Glover1
1Universität Heidelberg, Zentrum für Astronomie, Institut für Theoretische
Astrophysik, Albert-Ueberle-Str. 2, 69120 Heidelberg, Germany
2Department of Astronomy and Astrophysics, University of California, 1156 High
Street, Santa Cruz, CA 95064, USA
3Kavli Institute for Particle Astrophysics and Cosmology, Stanford University,
SLAC, Menlo Park, CA 94025, USA

###### Abstract

We employ the $\Delta$-variance analysis and study the turbulent gas dynamics
of simulated molecular clouds (MCs). Our models account for a simplified
treatment of time-dependent chemistry and the non-isothermal nature of the
gas. We investigate simulations using three different initial mean number
densities of $n_{0}=30,100$ and $300\>{\rm cm^{-3}}$ that span the range of
values typical for MCs in the solar neighbourhood. Furthermore, we model the
CO line emission in a post-processing step using a radiative transfer code. We
evaluate $\Delta$-variance spectra for centroid velocity (CV) maps as well as
for integrated intensity and column density maps for various chemical
components: the total, H2 and 12CO number density and the integrated intensity
of both the 12CO and 13CO ($J=1\rightarrow 0$) lines. The spectral slopes of
the $\Delta$-variance computed on the CV maps for the total and H2 number
density are significantly steeper compared to the different CO tracers. We
find slopes for the linewidth-size relation ranging from 0.4 to 0.7 for the
total and H2 density models, while the slopes for the various CO tracers range
from 0.2 to 0.4 and underestimate the values for the total and H2 density by a
factor of $1.5-3.0$. We demonstrate that optical depth effects can
significantly alter the $\Delta$-variance spectra. Furthermore, we report a
critical density threshold of $\sim 100\,$cm-3 at which the $\Delta$-variance
slopes of the various CO tracers change sign. We thus conclude that carbon
monoxide traces the total cloud structure well only if the average cloud
density lies above this limit.

###### keywords:

galaxies: ISM – ISM: clouds – ISM: molecules

##  1 Introduction

The interstellar medium (ISM) is dominated by highly turbulent motions, which
contribute to regulating stellar birth in molecular clouds (MCs) (see, e.g.
Mac Low & Klessen, 2004; Scalo & Elmegreen, 2004; Elmegreen & Scalo, 2004;
McKee & Ostriker, 2007; Hennebelle & Falgarone, 2012). On large scales,
turbulent gas motions support MCs against gravitational collapse. However, the
shocks associated with supersonic turbulence will create overdense regions on
small scales, which in turn may collapse and form stars and clusters. In
addition, ISM turbulence influences other physical processes such as the
chemical makeup of the gas and the associated heating and cooling processes,
the efficiency with which the external radiation field will be able to
penetrate into the interior of dense clouds, as well as the overall magnetic
field. Better understanding the role of turbulence in the ISM, therefore, is
of pivotal importance for many fields of modern astrophysics (see also the
lecture notes by Klessen & Glover, 2014).

Comparing theoretical simulations to observational measurements is very
difficult, since observations are always a complex convolution of the density
and the velocity field, affected by several other important physical processes
(e.g. by magnetic fields, the chemistry, stellar feedback, etc.). Furthermore,
it is often difficult to infer reliable physical parameters from observational
measurements. For example, observers rely on abundant tracers, e.g. 12CO, in
order to measure the amount of H2 gas in MCs. An obstacle with the line
emission of 12CO is that it becomes optically thick in dense cloud regions and
thus is not a good tracer for those regions anymore. To avoid this problem,
observers also employ other tracers to study the mass distribution within a
cloud, e.g. dust or 13CO, which is often optically thin. In addition,
molecular tracers like CS, HCN, HCO+ or NH3 can be used to study the gas mass
in very dense cloud regions. However, comparisons with simulations are needed
in order to study the influence of chemical inhomogeneities and optical depth
effects on observational measurements.

Lis et al. (1996) introduced maps of centroid velocities (CV) and showed that
those are to a certain degree sensitive to the underlying physics.
Accordingly, it is possible to use two-point statistics in order to recover
important information (e.g. optical depth effects, density and velocity
fluctuations, etc.) of an astrophysical system, e.g. by using centroid
velocity increment structure functions (see, e.g. Hily-Blant et al., 2008;
Federrath et al., 2010), spectral correlation functions (Rosolowsky et al.,
1999), the velocity channel analysis (Lazarian & Pogosyan, 2000, 2004), the
Principal Component Analysis (Heyer & Schloerb, 1997; Brunt & Heyer, 2002a) or
the $\Delta$-variance (see, e.g. Stutzki et al., 1998; Mac Low & Ossenkopf,
2000; Ossenkopf et al., 2001; Bensch et al., 2001; Ossenkopf et al., 2008a,
b). The latter provides a wavelet-based measure for characterizing structures
in astronomical datasets. For example, Stutzki et al. (1998) used the
$\Delta$-variance in order to measure the scaling behavior of structures in
observed images by analyzing the line emission of 12CO and 13CO, while Bensch
et al. (2001) studied the influence of white noise and beam smoothing on the
$\Delta$-variance spectra. Moreover, Ossenkopf et al. (2008b) applied the
$\Delta$-variance to interstellar turbulence and observations of MCs and
tested the capabilities of the method in a practical use by applying different
filter functions with different diameter ratios.

Several studies tried to reveal the influence of projection effects on
different statistical quantities, finding that various physical processes
might influence the projection of three-dimensional data on a two-dimensional
map on the sky. For example, Lazarian & Pogosyan (2004), Burkhart et al.
(2013b), Burkhart et al. (2014) and Bertram et al. (2015a) have shown that
optical depth effects might significantly alter the statistics of centroid
velocities due to different opacities of the gas tracers. Furthermore, the
turbulent driving as well as density and temperature fluctuations along the
line-of-sight also have a significant impact on the CV statistics (Lazarian &
Esquivel, 2003; Ossenkopf et al., 2006; Esquivel et al., 2007; Hily-Blant et
al., 2008). Much progress has been made during the last years in this context,
although we are still missing a consistent picture of turbulence theory in the
ISM.

In this paper we study high-resolution, 3D and time-dependent chemistry models
of hydrodynamical numerical simulations of the turbulent ISM. We use the
$\Delta$-variance method (see, e.g. Stutzki et al., 1998; Bensch et al., 2001;
Ossenkopf et al., 2008a) in order to study the structure of MCs and analyze
how different chemical tracers affect the $\Delta$-variance spectra. We also
vary the initial number density in order to explore how different densities
affect the statistics. Moreover, we perform a radiative transfer post-
processing and produce synthetic maps as well as position-position-velocity
(PPV) cubes of the 12CO and 13CO emission. We then compute maps of centroid
velocities, column densities and integrated intensities in order to analyze
how the chemical inhomogeneity and the variable opacity affect the various
$\Delta$-variance spectra.

In Section 2 we present our numerical simulations, introduce the radiative
transfer post-processing and explain the statistical methods. In Section 3 we
present our results, which are discussed in Section 4. Finally, we summarize
our findings and present our conclusions in Section 5.

##  2 Methods and simulations

The simulation data presented in the following sections are also used in
Bertram et al. (2014) and Bertram et al. (2015a). While Bertram et al. (2014)
use the technique of Principal Component Analysis (PCA) in order to
statistically analyze the turbulent flows in the spectral data cubes, Bertram
et al. (2015a) evaluate the slopes of the structure functions of 2D projected
centroid velocities and compute Fourier spectra. Here, we summarize the most
important aspects of our hydrodynamical simulations and the radiative transfer
post-processing.

###  2.1 Computational method

The simulations in this paper were performed using a modified version of the
Zeus-Mp MHD code (Norman, 2000; Hayes et al., 2006). We have embedded a
detailed atomic and molecular cooling function, described in Glover et al.
(2010) and Glover & Clark (2012), together with a simplified treatment of the
molecular gas chemistry. The chemical network 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. The network tracks the abundances of 9 species and follows 30
chemical reactions. We adopt the standard solar abundances of hydrogen and
helium. The abundances of carbon and oxygen are taken from Sembach et al.
(2000). We use $x_{\text{C}}=1.4\times 10^{-4}$ and $x_{\text{O}}=3.2\times
10^{-4}$, where $x_{\text{C}}$ and $x_{\text{O}}$ are the fractional
abundances by number of carbon and oxygen relative to hydrogen (for further
reading see also Glover et al., 2010; Glover & Mac Low, 2011). When we start
our simulations, the carbon is assumed to be singly ionised as C+, while the
hydrogen, helium as well as the oxygen are in atomic form. For further
information about the chemical model we refer the reader to Glover & Clark
(2012).

We run simulations with three different initial number densities $n_{0}$ of
the hydrogen nuclei, $n_{0}=300,100$ and $30\>{\rm cm^{-3}}$. The temperature
of the gas at the beginning of our runs is set to a constant value of 60 K. We
use an uniform weak magnetic field strength $B_{0}=5.85\mbox{ }\mu{\rm G}$,
which is initially oriented parallel to the $z$-axis of the computational
domain. We do not include self-gravity. The gas is initially uniform and
embedded in a periodic box with a side length of 20 pc. The turbulent
simulations are uniformly driven between wavenumbers $k=1$ and $k=2$ (Mac Low
et al., 1998; Mac Low, 1999) with a 3D rms velocity of $v_{\rm rms}=5\>{\rm
km\>s^{-1}}$. This value remains approximately constant throughout the whole
simulation period. For the dust-to-gas ratio, we adopt the standard local
value of 1:100 (Glover et al., 2010), and assume that the dust properties do
not vary with the gas density. The cosmic ray ionization rate was set to
$\zeta=10^{-17}\>{\rm s^{-1}}$. For the incident ultraviolet radiation field,
we adopt the standard parameterization of Draine (1978). This field has a
strength $G_{0}=1.7$ in Habing (1968) units, corresponding to an integrated
flux of $2.7\times 10^{-3}\,$erg cm-2s-1. We use a numerical resolution of
$512^{3}$ grid cells and discuss the influence of our limited numerical
resolution in more detail in Appendix A. Using this value for the resolution
as well as the scale of the total box, we can compute the cubic cell size and
obtain $\Delta x\approx 0.04\,$pc for our numerical simulations.

###  2.2 Radiative transfer post-processing

We use the radiative transfer code Radmc-3d111www.ita.uni-
heidelberg.de/$\sim$dullemond/software/radmc-3d/ (Dullemond, 2012) in order to
model the CO ($J=1\rightarrow 0$) line for both 12CO and 13CO. Furthermore, we
apply the Large Velocity Gradient (LVG) approximation (Sobolev, 1957) to
compute the level populations, as explained by Shetty et al. (2011a). Beyond
the line emission, Radmc-3d also accounts for the dust continuum emission,
which we subtract off before we analyze the statistical properties of the
synthetic maps. We use a number of 512 channels in velocity space for our
radiative transfer post-processing, corresponding to a spectral resolution of
$\sim 0.07\,$km s-1, $\sim 0.06\,$km s-1 and $\sim 0.05\,$km s-1 for the n300,
n100 and n30 model, respectively.

Our simulations do not explicitly track the abundance of 13CO (which would be
costly), and so we need a procedure to relate the 13CO number density to that
of 12CO. A common assumption is that the ratio of 12CO to 13CO is identical to
the elemental abundance ratio of 12C to 13C (see, e.g. Roman-Duval et al.,
2010). In the majority of our analysis, we make the same assumption and set
the 12CO to 13CO ratio to a constant value, $R_{12/13}=50$, which we use to
generate a first set of 13CO column density maps. However, the effects of
chemical fractionation (Watson et al., 1976) and selective photodissociation
of 13CO (see, e.g. Visser et al., 2009) can significantly alter the value of
$R_{12/13}$ within the cloud (Röllig & Ossenkopf, 2013; Szűcs et al., 2014).
Therefore, we also explore the effect that this may have on the statistics of
the $\Delta$-variance, using numerical results of Szűcs et al. (2014). Szűcs
et al. (2014) give a numerical fitting formula, which relates the ratio
$R_{12/13}$ to the 12CO column density in the cloud. Thus, we also produce a
second set of maps using their routine and compute variable 13CO number
densities by dividing the 12CO number densities by the mean value of
$R_{12/13}$ from the fitting formula for each line-of-sight. Finally, we can
use these 13CO number densities in order to compute the 13CO emission in the
same fashion as in our constant ratio models. The differences in the
$\Delta$-variance spectra of the two sets of 13CO 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 Section 3.3 below.

The radiative transfer calculation yields position-position-velocity (PPV)
cubes of brightness temperatures $T_{B}$, which are related to the intensity
via the Rayleigh-Jeans approximation,

$T_{B}(\nu)=\left(\frac{c}{\nu}\right)^{2}\frac{I_{\nu}}{2k_{B}},$ (1)

where $I_{\nu}$ is the specific intensity at frequency $\nu$ and $k_{B}$ the
Boltzmann constant. We will refer to these as the “intensity” models. In
analogy, we construct centroid velocity and column density maps out of the PPP
simulation data of the density and velocity field, as described in the
following Section 2.3. We will refer to these as the “density” models,
following the notation already used in Bertram et al. (2014) and Bertram et
al. (2015a).

###  2.3 Centroid velocity, integrated intensity and column density maps

One part of our study is based on the statistics of centroid velocities. The
centroids $C(x,y)=C(\textbf{r})$ are defined as

$C(\textbf{r})=\frac{\int F(\textbf{r},z)v_{z}(\textbf{r},z)\text{d}z}{\int
F(\textbf{r},z)\text{d}z},$ (2)

where the variable $v_{z}(\textbf{r},z)$ is the velocity component along the
line-of-sight (e.g. the $z$-direction). Thus, $C(\textbf{r})$ is a map of
line-of-sight projected velocities. The quantity $F(\textbf{r},z)$ is a
statistical weight. It denotes either the underlying density field or the
brightness temperatures from the PPV cubes. As shown by Lis et al. (1996),
centroid velocity maps are very sensitive to the underlying physics of a MC,
which means that they can be used in order to infer important information
about the astrophysical system. Furthermore, CV maps can be easily computed
from observational data that rely on spectral measurements. These can be then
translated into maps of centroid velocities via equation (2).

Beside the CV maps, we also consider maps of integrated intensity and column
density. We calculate the former along a given line-of-sight (e.g. the
$z$-direction) via

$W(\textbf{r})=\int T_{B}(\textbf{r},z)\,\text{d}v,$ (3)

where $T_{B}$ is the brightness temperature, as introduced in Section 2.2. The
latter is given by

$N(\textbf{r})=\int n(\textbf{r},z)\,\text{d}z,$ (4)

where $n$ is the number density of an individual chemical component in a given
cell in the simulation box.

###  2.4 The $\Delta$-variance method

Figure 1: Example of a $\Delta$-variance spectrum, plotted as a function of
spatial scale, normalized by the total box size and averaged over all
available time snapshots. The inferred slope $\alpha$ is also indicated (solid
line). The $\Delta$-variance is computed on a CV map, where the velocities are
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-$\sigma$ spatial and temporal variations.

We use the $\Delta$-variance tool developed by Ossenkopf et al.
(2008a)222http://hera.ph1.uni-koeln.de/ ossk/Myself/deltavariance.html. The
$\Delta$-variance method measures the variance in a structure $S(\textbf{r})$
(in our case the maps of centroid velocities, integrated intensities and
column densities) on a given spatial scale $\ell$, by filtering the dataset
with a spherically symmetric up-down-function of size $\ell$. It is given by

$\sigma_{\Delta}^{2}(\ell)=\Bigl{\langle}\Big{(}S(\textbf{r})\ast\bigodot_{\ell}(\textbf{r})\Big{)}^{2}\Bigl{\rangle}_{\textbf{r}},$
(5)

where the average is computed over all positions $\textbf{r}=(x,y)$ on the
sky. The symbol $\ast$ stands for a convolution and $\bigodot_{\ell}$
describes the filter function. In this paper, we use a Mexican hat with a
diameter ratio of 1.5. However, we have also analyzed the impact of the filter
function and the diameter ratio on the slope values of the spectra, e.g. by
using a French hat with a diameter ratio of 3.0. We find that our results do
not significantly depend on the specific choice of the filter function and the
diameter ratio (see Appendix B). The differences between the individual filter
functions are described in more detail in Ossenkopf et al. (2008a).

In this study we compute the $\Delta$-variance for each possible line-of-sight
direction $x$, $y$ and $z$. As shown by Esquivel & Lazarian (2005) and
Burkhart et al. (2014), the statistics of velocity centroids are very
sensitive to the direction of the magnetic field in the regime of sub-Alfvénic
Mach numbers. Although we use a weak magnetic field in the $z$-direction, the
turbulence in our simulations is trans-Alfvénic or mildly super-Alfvénic (see
Table 1). The field lines are essentially dragged along with the turbulent
flow, with the result that the turbulence remains approximately isotropic.
Hence, we do not find significant variations of the $\Delta$-variances along
the different directions. We therefore average all $\Delta$-variances of the
three line-of-sights. Finally, power-laws of the form

$\sigma_{\Delta}^{2}(\ell)\propto\ell^{\alpha}$ (6)

were fit to the $\Delta$-variance spectra, where $\alpha$ denotes the slope of
the power-law. To calculate the scaling exponents, we use a fitting range from
$1/10$ to $1/4$ of the total box size ($0.1\lesssim\ell/D\lesssim 0.25$), as
constrained by Federrath et al. (2010) and Konstandin et al. (2012). For a box
with $D=512$ grid cells for each side this translates to 51 and 128 cells,
corresponding to a physical scale of $\sim 2-5\,$pc in the simulation domain.
Extending the fitting range to scales above this limit is complicated since
the simulations are driven on large scales, which would significantly bias our
results. Additionally, the scales below our limit are influenced by the
numerical resolution and the bottleneck effect (Kritsuk et al., 2007;
Konstandin et al., 2015), which is an accumulation of kinetic energy caused by
the viscosity of the fluid before it dissipates into thermal energy.
Furthermore, as demonstrated by Konstandin et al. (2015), we also note that a
clear and unique identification of the fitting range in numerical simulations
with our resolution is complicated. Thus, the specific choice of the fitting
range introduces some degree of uncertainty in the $\Delta$-variance slopes,
which we conservatively estimate to be of the order of $\pm 0.1$.

Fig. 1 shows an example of a $\Delta$-variance spectrum computed using maps of
centroid velocities inferred from our n300 run. The velocities are weighted by
the total density field and the spectrum is averaged over all available time
snapshots. We also show the fitting slope $\alpha$ for this spectrum. This
slope can be translated into a linewidth-size relation, given by
$\sigma_{\Delta}(\ell)\propto\ell^{\gamma}$, with a scaling exponent
$\gamma=\alpha/2$. The latter equation relates the velocity fluctuations to
their spatial scale, which is often referred to as Larson’s law (Larson,
1981). This relation suggests that turbulence plays an important role in the
process of star formation and predicts a turbulent energy cascade within the
cloud, as proposed by Kolmogorov (1941). Hence, the slopes $\gamma$
characterize the turbulent velocity hierarchy in our clouds and can be
directly compared to values derived from spectral observations of MCs,
typically ranging from $\gamma\approx 0.2-0.9$ (see, e.g. Larson, 1981;
Solomon et al., 1987; Brunt & Heyer, 2002b; Heyer & Brunt, 2004; Heyer et al.,
2006; Sun et al., 2006; Hily-Blant et al., 2008; Roman-Duval et al., 2011;
Schneider et al., 2011; Russeil et al., 2013; Elia et al., 2014; Alves de
Oliveira et al., 2014; Johnston et al., 2014). For further reading about the
theory of turbulence in astrophysics and fluids, we refer the reader to
Burgers (1948), Benzi et al. (1993), Frisch & Kolmogorov (1995), She & Leveque
(1994), Dubrulle (1994), Boldyrev et al. (2002) or Mac Low & Klessen (2004).

Furthermore, the scaling exponents $\alpha$ of the $\Delta$-variance spectra
are related to the power-law exponents $\beta$ of the corresponding Fourier
spectra via $\alpha=\beta-1$ for 3D data (Stutzki et al., 1998), where the
power-law spectrum in Fourier space is defined as

$P(k)\propto k^{-\beta}.$ (7)

In this equation, $k=2\pi/\ell$ denotes the wavevector. In analogy, for a
power spectrum of a 2D image, the $\Delta$-variance is related to the power
spectrum via $\alpha=\beta-2$.

We also note that there are many other techniques to measure structural
density and velocity fluctuations as a function of spatial scale, e.g. by
computing structure functions (see, e.g. Hily-Blant et al., 2008; Federrath et
al., 2010) or the spectral correlation function (Rosolowsky et al., 1999), by
using the velocity channel analysis (Lazarian & Pogosyan, 2000, 2004) or the
Principal Component Analysis (Heyer & Schloerb, 1997; Brunt & Heyer, 2002a;
Bertram et al., 2014).

##  3 Results

Model name | n300 | n100 | n30
---|---|---|---
Mean density $[$cm${}^{-3}]$ | 300 | 100 | 30
Resolution | 5123 | 5123 | 5123
Box size [pc] | 20 | 20 | 20
$t_{\text{end}}$ [Myr] | 5.7 | 5.7 | 5.7
$\langle\mathcal{M}_{\text{s}}\rangle$ | 10.6 | 6.8 | 5.1
$\langle\mathcal{M}_{\text{A}}\rangle$ | 1.5 | 1.1 | 1.0
$\sigma_{\rho}/\langle\rho\rangle$ | 3.0 | 3.0 | 4.6
$\langle x_{\text{H}_{2}}\rangle_{\text{mass}}$ | 0.98 | 0.78 | 0.61
$\langle n_{\text{H}_{2}}\rangle_{\text{vol}}$ $[$cm${}^{-3}]$ | 140 | 37 | 9
$\langle n_{\text{H}_{2}}\rangle_{\text{mass}}$ $[$cm${}^{-3}]$ | 1456 | 447 | 264
$\langle n_{\text{CO}}\rangle_{\text{vol}}$ $[$cm${}^{-3}]$ | $2.2\times 10^{-2}$ | $1.7\times 10^{-3}$ | $1.4\times 10^{-4}$
$\langle n_{\text{CO}}\rangle_{\text{mass}}$ $[$cm${}^{-3}]$ | 0.3 | $4.3\times 10^{-2}$ | $2.3\times 10^{-2}$
$\langle T\rangle_{\text{vol}}$ $[$K$]$ | 35 | 68 | 223
$\langle T\rangle_{\text{mass}}$ $[$K$]$ | 13 | 26 | 57
$\langle N_{\text{H}_{2}}\rangle$ $[$cm${}^{-2}]$ | $8.6\times 10^{21}$ | $2.3\times 10^{21}$ | $5.3\times 10^{20}$
$\langle N_{\text{CO}}\rangle$ $[$cm${}^{-2}]$ | $1.3\times 10^{18}$ | $1.1\times 10^{17}$ | $8.7\times 10^{15}$
Table 1: Overview of our different models with some characteristic values
measured for the last time snapshot. From top to bottom we list: mean number
density, resolution, box size, time of the last snapshot, mean sonic Mach
number, mean Alfvén Mach number, ratio of density dispersion and mean density,
mass-weighted mean abundances of H2 (i.e. the percentage of atomic hydrogen
that has been converted to H2, see Glover & Mac Low, 2011), mean volume- and
mass-weighted H2 and CO number densities, mean volume- and mass-weighted
temperature and mean column densities of H2 and CO.

Figure 2: Logarithmic maps of column density as well as velocity-integrated
intensity along the $z$-direction for our n30 (left column), n100 (middle
column) and n300 models (right column). From top to bottom we show the
different chemical components: total, H2 and CO column density as well as the
integrated intensity of 12CO and 13CO in the $J=1\rightarrow 0$ transition.
Each side of the simulation domain has a length of 20 pc. Note that the
velocity field of the n30 model uses a different turbulent random seed than
the n100 and the n300 model. Furthermore, we caution the reader that our color
bars use a different scaling in the individual plots.

Figure 3: Same as Fig. 2, but with maps of centroid velocities (CV). Black
areas in the 12CO and 13CO intensity map of the n30 model denote regions where
the brightness temperatures are zero along the line-of-sight. We mask these
regions, because no proper centroid velocities can be computed there.

We perform numerical simulations and apply radiative transfer post-processing
to our data in order to analyze the influence of chemical inhomogeneities and
optical depth effects on the $\Delta$-variance analysis. Table 1 gives an
overview of our numerical models. Mass- and volume-weighted quantities are
defined via $\langle f\rangle_{\text{mass}}=\sum
f\rho\text{d}V/\sum\rho\text{d}V$ and $\langle f\rangle_{\text{vol}}=\sum
f\text{d}V/\sum\text{d}V$, respectively. For more information about the H2 and
CO distributions produced in this kind of turbulent simulation as well as
integrated intensity and column density PDFs, we refer the reader to Glover et
al. (2010) and Shetty et al. (2011a).

Fig. 2 and 3 show the column density and integrated intensity maps as well as
the centroid velocity maps computed via equations (2), (3) and (4) for all
models and chemical components. Regarding Fig. 2, we find that the total and
H2 density models show similar structures on all spatial scales, which is
because most of the hydrogen is in molecular form at this time (see also Table
1). The CO column densities also trace similar structures, but span a much
wider range of values, demonstrating that carbon monoxide has very low
abundances along the low column density LoS. Furthermore, the 12CO and 13CO
intensity maps also largely reflect the distribution of the carbon monoxide
gas. We find that the intensity maps are much smoother and span a smaller
range of values ($4-7$ orders of magnitude in integrated intensity compared to
$7-18$ orders of magnitude in column density), which is due to the fact that
the $J=1\rightarrow 0$ line of 12CO is easily excited and can be bight even in
low-density cloud regions. In the case of the 13CO intensity maps, the
abundance is lower by a factor of $R_{12/13}=50$ compared to the 12CO
intensity models. Nevertheless, 13CO is optically thin and so the peaks in the
13CO maps coincide with those of the CO density models.

###  3.1 $\Delta$-variance analysis of the CV maps

Figure 4: $\Delta$-variance spectra as a function of the spatial scale, averaged over all available time snapshots and the different lines-of-sights $x$, $y$ and $z$. The spatial scale is normalized by the total box size. From top to bottom: spectra for our three different density models, i.e. n300, n100 and n30, computed for the centroid velocity (CV) maps. In each plot we show the $\Delta$-variance spectra for our various chemical components, i.e. for the total density, H2 and CO density as well as for the 12CO and 13CO intensity. We use a fitting range from $1/10$ to $1/4$ of the total box size ($0.1\lesssim\ell/D\lesssim 0.25$), as constrained by Federrath et al. (2010) and Konstandin et al. (2012), which is indicated by a horizontal dashed line and is the same in each density model. The different power-law functions in the fitting range are indicated with a black solid line on each spectrum. Furthermore, we list the slopes $\alpha$ and their errors from a $\chi^{2}$-fit in each plot. | Total density | H2 density | CO density |  12CO intensity |  13CO intensity
---|---|---|---|---|---
| $\alpha$ | $\alpha$ | $\alpha$ | $\alpha$ | $\alpha$
n300 | $1.10\pm 0.02$ | $1.06\pm 0.02$ | $0.66\pm 0.03$ | $0.75\pm 0.03$ | $0.73\pm 0.02$
n100 | $0.83\pm 0.01$ | $0.70\pm 0.01$ | $0.28\pm 0.01$ | $0.30\pm 0.01$ | $0.31\pm 0.02$
n30 | $1.29\pm 0.01$ | $1.27\pm 0.03$ | $0.86\pm 0.02$ | $0.85\pm 0.02$ | $0.84\pm 0.02$
Table 2: Spectral slope of the $\Delta$-variance spectrum, computed using
centroid velocities weighted by the indicated quantity (see also Fig. 4). The
errors are computed by the $\chi^{2}$ fitting method. The slope $\alpha$ is
related to the $\Delta$-variance via
$\sigma_{\Delta}^{2}(\ell)\propto\ell^{\alpha}$ in a linear regime. The values
$\alpha$ can also be used to compute spectral slopes $\gamma=\alpha/2$ for a
linewidth-size relation $\sigma_{\Delta}(\ell)\propto\ell^{\gamma}$, readily
comparable to observational measurements.

We compute the $\Delta$-variance for all models (n300, n100 and n30) and all
chemical tracers: the total density, H2 density, CO density, 12CO and 13CO
intensity. We average all spectra from 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. Fig. 4 shows
$\Delta$-variance spectra as a function of spatial scale for all density
models and chemical tracers, computed for the centroid velocity (CV) maps, as
introduced in Section 2.3. Table 2 summarizes all slope values $\alpha$
obtained from a $\chi^{2}$-fit for the $\Delta$-variance of the velocity
field, $\sigma_{\Delta}^{2}(\ell)\propto\ell^{\alpha}$, which are also listed
in each plot in Fig. 4.

We generally find significant differences between the various models and
tracers. In our n300 model, we obtain similar slopes for the total density and
the H2 density (see Table 2). Regarding the spectra of those two cases, we
find an excellent agreement within the values of $\sigma_{\Delta}^{2}$ for
both models, as shown in the top plot of Fig. 4. For the low-density model n30
(bottom plot of Fig. 4), we obtain a large discrepancy between the
$\sigma_{\Delta}^{2}$ values of total density and the H2 density model. This
is because in this simulation the fraction of molecular gas is much smaller
than in the higher density runs and so H2 no longer follows the total gas
density (Glover et al., 2010). The H2 density is therefore more inhomogeneous
than the gas density, and as a result the centroid velocities weighted by H2
fluctuate more on all scales than those weighted by the total density. The
n100 model makes up an intermediate case between the n300 and the n30 density
models. In this model, the correlation between the total and the H2 density is
worse than in the n300 case, but still better than in the n30 model.

Comparing the total and H2 density models to the various CO tracers, we find
that the slopes of the former are significantly steeper compared to the slopes
of the latter, independent of the density (see Table 2). However, the
$\Delta$-variance spectra of the various CO tracers show a good agreement with
each other over nearly all spatial scales. Furthermore, we always find the
$\sigma_{\Delta}^{2}$ values of the CO tracers to be significantly larger than
those of the total and H2 density. This is because CO is mainly located in
denser regions of the cloud (see, e.g. Shetty et al., 2011a,b), leading to
higher density contrasts compared to the total and H2 density and thus to
larger variances in velocity space.

###  3.2 $\Delta$-variance analysis of the intensity and column density maps

Figure 5: Same as Fig. 4, but computed using maps of integrated intensity and column density. In this plot, we normalize each curve by the corresponding value $\sigma_{0}^{2}$ measured at an arbitrary spatial scale of 5% of the total box, in order to better compare the spectra computed for the intensity and column density maps with each other. For a better visualization, we also show a horizontal dashed line at a value of $\sigma_{\Delta}^{2}/\sigma_{0}^{2}=1$ in each plot. Furthermore, for clarity, we do not plot the error bars of the individual spectra here. The slopes $\alpha$ and their errors from a $\chi^{2}$-fit are listed in each plot for all our chemical components. | Total density | H2 density | CO density |  12CO intensity |  13CO intensity
---|---|---|---|---|---
| $\alpha$ | $\alpha$ | $\alpha$ | $\alpha$ | $\alpha$
n300 | $0.34\pm 0.02$ | $0.32\pm 0.02$ | $0.37\pm 0.02$ | $0.23\pm 0.02$ | $0.22\pm 0.01$
n100 | $0.20\pm 0.01$ | $0.10\pm 0.01$ | $-0.18\pm 0.01$ | $-0.04\pm 0.01$ | $-0.06\pm 0.01$
n30 | $0.51\pm 0.01$ | $0.39\pm 0.03$ | $-0.98\pm 0.04$ | $-0.35\pm 0.04$ | $-0.72\pm 0.01$
Table 3: Same as Table 2, but for the $\Delta$-variance computed directly from
maps of integrated intensity and column density, as defined in equation (3)
and (4). The corresponding $\Delta$-variance spectra are shown in Fig. 5.

Beyond the $\Delta$-variance analysis of centroid velocities, we can also
apply the same method to any other quantities defined on the $x$-$y$ plane. In
particular, it is possible to carry out a similar analysis for the total, H2
and CO column densities and the 12CO and 13CO integrated intensities. The
results of this analysis are illustrated in Fig. 5. In contrast to Fig. 4, we
normalize each curve by the corresponding $\Delta$-variance value
$\sigma_{0}^{2}$ measured at an arbitrary spatial scale of 5% of the total
box, in order to better compare the spectra computed for the intensity and
column density maps with each other. Furthermore, we also show a horizontal
dashed line at a value of $\sigma_{\Delta}^{2}/\sigma_{0}^{2}=1$ in each plot.
Table 3 summarizes all slope values $\alpha$ obtained from a $\chi^{2}$-fit
for the individual $\Delta$-variance spectra.

As for the CV maps, we find significant differences between the spectra of the
various models and tracers computed using the maps of integrated intensity and
column density directly. In the high-density n300 model, we find a good
correlation between the total density and the H2 density (see also the slopes
in Table 3), for the same reasons as described in Section 3.1. Similarly, the
low-density n30 model shows again the largest discrepancy. This is because the
fraction of molecular gas is significantly smaller than in the n300
simulation, leading to different values of $\alpha$. The n100 model again
defines an intermediate case between the n300 and the n30 density runs. In
this model, the H2 density better correlates with the total density than in
the low-density n30 run, but still worse than in the high-density n300 run.

Regarding the simulations of the various CO tracers, we obtain very different
scaling properties compared to the total and the H2 density simulations,
depending on the mean density. In the high-density n300 model, we find a good
correlation between the two CO intensity cases, while the slopes of the total
and CO density models are significantly steeper. However, both the total and
H2 density as well as the various CO tracer models show the same trend in
$\sigma_{\Delta}^{2}$ over all spatial scales, having $\alpha>0$. This
situation changes in the low-density n30 model, and we find strong differences
between the CO tracers and the total and H2 density runs over all spatial
scales. While the structure in the computational domain on larger scales for
the total and H2 column density can still be described by a slope $\alpha>0$,
the $\Delta$-variance spectra for the CO tracers have negative slopes in this
density run (see Table 3). A similar situation as in the low-density n30
simulation can be seen in the intermediate n100 density model. While we again
measure $\alpha>0$ for the total and H2 density cases, we obtain
$\alpha\lesssim 0$ for the CO tracers, meaning that the spatial distribution
of observable carbon monoxide is completely different compared to the
distribution of the H2 gas in the n30 and the n100 density simulations.

Looking at Fig. 5, we see that the reason that we obtain negative values for
$\alpha$ for the CO tracers in the low-density runs is that the corresponding
$\Delta$-variance spectra turn over at relatively low spatial scales. This
indicates that the CO in these simulations has a characteristic spatial scale
of approximately 2% of the total box size, corresponding to around 0.4 pc.
This is consistent with what we see in the CO column density projections and
integrated intensity maps shown in Fig. 2. In the low-density n30 simulation,
CO is sufficiently well shielded to resist photodissociation only in the dense
filaments formed by turbulent compressions. These structures are thin and
contribute strongly to the $\Delta$-variance on scales comparable to their
width. Moreover, the fact that this characteristic scale is roughly the same
in all CO spectra demonstrates that the radiative transfer post-processing
does not significantly change the characteristic spatial cloud signatures in
the $\Delta$-variance analysis. Furthermore, 12CO is optically thick and thus
we always find slightly steeper $\alpha$ slopes for the 12CO intensity tracers
compared to the optically thin 13CO intensity tracers.

###  3.3 Variation of the abundance of 13CO

We also analyze the variation of the abundance of 13CO, because our
simulations only follow the chemistry of 12CO, but not of its isotope 13CO.
Therefore, we produce a set of 13CO emission maps, using a spatially varying
$R_{12/13}$ generated following the prescription given in Section 2.2 and in
Szűcs et al. (2014). After applying the radiative transfer post-processing to
these maps, we compute the maps of the CV and the integrated intensity and
compare the $\Delta$-variance spectra with each other. The results agree with
those as described in Bertram et al. (2015a) for the centroid velocity
increment structure functions. Although we encounter slight variations on
small scales for both the CV and the integrated intensity maps, the spectra of
the $\Delta$-variance do not change significantly as we vary $R_{12/13}$.
Thus, the use of a constant 12CO/13CO ratio is sufficient for obtaining proper
$\Delta$-variance spectra for 13CO.

##  4 Discussion

###  4.1 Interpreting the $\Delta$-variance spectra computed using the CV maps

Regarding the spectra computed on the CV maps in Fig. 4, we find a good
correlation between the total density and the H2 density in the n300 model.
This is because $\sim 98$% of the initial atomic hydrogen is in molecular form
at this time. Thus, the spectra of the $\Delta$-variance for the CV maps are
primarily dominated by the H2 mass in this case. Vice versa, if we analyze the
CV spectra of the n30 model, we find a larger discrepancy between the total
density and the H2 density, which is because only $\sim 61$% of the initial
atomic hydrogen is in molecular form at this time. The n100 model is an
intermediate case, with this value being $\sim 78$%. Furthermore, we obtain a
good correlation between the various CO tracer components. This indicates,
that the turbulence statistics are similar for all three cases and that the
impact of the radiative transfer post-processing on the $\Delta$-variance
analysis is limited. This is similar to the conclusion presented in Bertram et
al. (2015a) based on the analysis of CV increment structure functions.

In general, we find that the slopes $\alpha$ for the total and H2 density
models are significantly steeper than the slopes for the different CO tracer
models (see the values in Table 2). This indicates that these components have
a higher relative amount of structures on larger scales compared to the
different CO tracers. Furthermore, as shown by Glover et al. (2010) and Shetty
et al. (2011a,b), CO is primarily a good tracer of dense and very compact
regions in a cloud. Thus, this leads to less turbulent velocity structures on
larger scales and hence to flatter slopes compared to those of the total and
H2 density. However, we find that the different slopes $\alpha$ of the various
CO tracers underestimate the slopes of the total and H2 density by a factor of
$\sim 1.5-3.0$ (see Table 2).

Furthermore, we caution the reader that it is difficult to infer a clear
dependence of the slopes on the mean ISM density. It is likely that the
statistical measures we derive from our numerical simulations are also
sensitive to the specific realization of the turbulent velocity field. Since
we are studying flows which are driven on large scales, variance effects can
become important, and the statistical properties depend on the random
orientation of the turbulent modes as well (see, e.g. Klessen, 2000; Klessen
et al., 2000). Hence, in order to obtain slopes that properly converge with
the density, we speculate that a large number of simulations with various
turbulent realizations would be needed in order to calculate reliable average
values (see, e.g. the PCA analysis of the statistics using different
realizations of the turbulent velocity field in Bertram et al., 2014). This is
prohibited by the high computational costs of the individual simulations, and
hence we only focus on one example in this paper, which is enough to
illustrate basic trends of the $\Delta$-variance statistics.

We also note that a direct comparison of our 2D $\Delta$-variance statistics
to 3D turbulence models is difficult, since all our 2D maps are a complex
convolution of the 3D density field (or the brightness temperatures) with the
3D velocity field. Various different physical processes can influence the CV
statistics, as shown by previous studies. For example, Lazarian & Pogosyan
(2004), Burkhart et al. (2013a) and Burkhart et al. (2014) analyzed the impact
of the sonic Mach number on the CV statistics, finding that it can
significantly alter the results. Moreover, Lazarian & Esquivel (2003),
Ossenkopf et al. (2006), Esquivel et al. (2007), Hily-Blant et al. (2008) and
Federrath et al. (2010) studied the effects of the turbulent driving as well
as temperature and density fluctuations on the CV statistics, also finding
significant differences in the statistics and in the inferred CV slopes.
Nevertheless, although a direct comparison of the 2D statistical quantities to
3D measures is complicated, we can safely use the 2D $\Delta$-variance
analysis in order to work out statistical trends measured in the spectra of
the individual models.

###  4.2 Interpreting the $\Delta$-variance spectra computed on the intensity
and column density maps

In Section 3.1 and 3.2 we already justified the correlation between the total
and the H2 density runs in all our density models and compared the different
slopes $\alpha$ from Table 3 with each other. Furthermore, we established that
the observed CO gas distribution in the low-density n30 model does not reflect
the spatial distribution of molecular hydrogen very well. This is because CO
is mainly located in the dense filaments, while H2 is more space-filling and
distributed over the total cloud, owing to its greater ability to resist
photodissociation. Consequently, the $\Delta$-variance spectra of the CO
tracers peak at a scale corresponding to the width of these structures. Hence,
if we apply the $\Delta$-variance analysis to the maps of integrated intensity
or column density, we can obtain important information about the turbulently
created high-density peaks within the MC (see also the discussion in Ossenkopf
et al., 2001). However, we also note that the physical connection of our
slopes $\alpha$ from Table 3 to observational measurements is complicated,
since observations typically probe smaller spatial scales, which our
simulations are not sensitive to due to the limited numerical resolution.

If we compare the CO tracers with the total and H2 density in the n300 model,
we find that the carbon monoxide better reproduces the statistical trends of
molecular hydrogen in this high-density run compared to the low-density n30
model. In the n300 run, we generally obtain $\alpha>0$ for all chemical
components, i.e. we find more cloud structures on larger spatial scales. In
contrast to the low-density n30 model, this is because carbon monoxide is not
only confined to small dense filaments, but instead is distributed over the
whole molecular cloud. In particular, a significant amount of CO gas can also
be observed between the numerous dense cloud regions. Consequently, we see the
$\Delta$-variance spectra peaking at the largest scales, consistent with
observational efforts and previous work on this field (see, e.g. Stutzki et
al., 1998; Ossenkopf et al., 2001; Bensch et al., 2001; Sun et al., 2006;
Schneider et al., 2011; Russeil et al., 2013; Alves de Oliveira et al., 2014;
Elia et al., 2014).

However, regarding the general trend of the slopes $\alpha$ from the n30 to
the n300 via the n100 model, we find that the values of $\alpha$ for the CO
tracers change sign above a critical density of $\sim 100\,$cm-3. This is
about the number density at which all values $\alpha$ become positive. Thus,
we conclude that carbon monoxide traces our total cloud structure well only if
the average cloud density lies significantly above a critical threshold of
$\sim 100\,$cm-3. If the mean density in the cloud is significantly smaller
than this limit, the observable CO gas does not properly trace the statistical
properties of the H2 gas in the cloud. We speculate that one may also see a
similar switch from CO tracing all of the structures to only the dense cores
and filaments if one increases the incident field strength or decreases the
metallicity, as in both cases, this makes it much easier to photodissociate
the diffuse CO (see, e.g. Glover & Clark 2012; Bertram et al. 2015b,
submitted). On the contrary, if the mean density in the cloud is too high, CO
becomes optically thick and so we suspect that there should also exist an
upper density limit above which CO does not properly trace structures of the
cloud anymore. We leave such an analysis for further investigations.

###  4.3 Comparison of the $\Delta$-variance to other statistical tools
applied to our simulation data

We have previously applied several other statistical methods to the same set
of simulation data in order to study the structural behavior of MCs. For
example, Bertram et al. (2014) applied Principal Component Analysis (PCA) to
the same data used in this paper, while Bertram et al. (2015a) carried out a
similar study using centroid velocity increment structure functions (CVISF)
and Fourier spectra.

Comparing the results of the different methods with each other, we generally
find consistent results between the different statistical analysis methods.
For example, the slopes of the CVISF for the total density and H2 density
models are also significantly steeper compared to the various CO tracers,
which is in good agreement to the results of the $\Delta$-variance analysis in
this paper. The same holds for the relative scaling of the different CO
models, which is the same in the analysis of the CVISF and the
$\Delta$-variance. Interestingly, we also find a close relation between the
PCA structure analysis and the spectra of the $\Delta$-variance for the CO
tracers. The PCA method does not find any structures on larger scales for the
CO tracers in the n30 model as well as for some single CO tracers in the n100
model. This result can also be reproduced in the different CO spectra of the
$\Delta$-variance analysis, i.e. where the gradient of the individual spectra
become negative at a characteristic scale of $\sim 0.4\,$pc, as described in
Sections 3.2. This is the situation where small clumpy structures of CO gas
dominate the overall composition of carbon monoxide in the MC. In this case,
the $\Delta$-variance shows that the most dominant CO structures in the
simulation domain can be found at small scales, leading to completely missing
principal components on scales larger than $\sim 1\,$pc, as presented in Fig.
3 in Bertram et al. (2014). Hence, applying the $\Delta$-variance analysis to
maps of integrated intensities or column densities gives an idea about
characteristic spatial scales in the cloud of interest. The $\Delta$-variance
thus provides a good statistical tool in order to study the relative gas
distribution on various cloud scales.

These results indicate that all these various statistical methods are
connected to each other and that characteristic structural properties of the
MC should be traced by each of them individually. Hence, all various
statistical methods have proven to yield reasonable results for the structure
analysis of MCs. Nevertheless, the advantage of the $\Delta$-variance method
is that it is fast and easy to implement, while the computation of structure
functions and principal components is more expensive. Thus, the
$\Delta$-variance analysis provides a useful and adequate tool in order to
quickly study the internal structure of a cloud. However, the advantage of the
structure function analysis is that their results can be easily compared to
theoretical models of turbulence.

###  4.4 Previous studies about $\Delta$-variances computed using CV maps

Several studies tried to compute the $\Delta$-variance spectra of CV maps in
the past and estimated turbulent slope values $\alpha$ from observational
measurements. For example, Ossenkopf et al. (2008b) used optically thick 12CO
($J=1\rightarrow 0$) maps of the Polaris Flare to compute centroid velocity
maps and to estimate the slope $\alpha$. The Polaris flare is an archetype of
a cloud midway between the diffuse and the molecular phases (Heithausen &
Thaddeus, 1990; Meyerdierks & Heithausen, 1996; Falgarone et al., 1998; André
et al., 2010; Miville-Deschênes et al., 2010). It is supposed to have a low
average density and thus can be compared to our n30 low-density model.
Adopting the same 12CO tracer, we find a value of $\alpha=0.85\pm 0.02$, in
agreement with the estimate of $\alpha\approx 0.9$ found by Ossenkopf et al.
(2008b).

Moreover, Ossenkopf & Mac Low (2002) computed the slope $\alpha$ for
hydrodynamic supersonic simulations driven at different wavenumbers. They find
one power-law range for all models and obtain $\alpha\approx 1.0$, which fits
into our range of power-law slopes $\alpha$ from $0.8-1.3$ for the total
density models. Furthermore, Federrath et al. (2009) computed
$\Delta$-variance slopes for numerical simulations with both solenoidal and
compressive forcing, finding a similar range of $\alpha$ values from $0.8-1.4$
for the turbulent velocity field. However, these simulations only use
supersonic isothermal turbulence, while our runs include more complex physics,
e.g. a chemical network, heating and cooling, various initial number densities
or the coupling to the radiation field. Thus, if we compare the total
variations of slopes $\alpha$ in their and in our models, we find that the
influence of the different forcing methods on the slopes $\alpha$ in their
simulations is large, while the impact of our additional physical effects
(varying density, optical depth effects, etc.) on the slope values remains
comparatively small.

###  4.5 Previous studies about $\Delta$-variances computed using intensity
and column density maps

As well as computing $\Delta$-variance spectra for maps of centroid
velocities, several studies also applied the $\Delta$-variance to maps of
integrated intensity or column density. For example, Bensch et al. (2001)
computed $\Delta$-variance spectra on maps of velocity-integrated intensity
for 12CO and 13CO ($J=1\rightarrow 0$) for various MCs in the Galaxy. Stutzki
et al. (1998) analyzed the $\Delta$-variance of an observed 12CO
($J=1\rightarrow 0$) image of the Polaris flare as a whole, while Ossenkopf et
al. (1998) studied the intensity map of one of its subclouds, MCLD 123.5+24.9.

Comparing all the various $\Delta$-variance spectra in Bensch et al. (2001) or
Stutzki et al. (1998) with each other, we find a value $\alpha>0$ for each of
them in the given fitting range, i.e. the spectra increase with increasing
spatial lag. At first sight, this is in contradiction to the results that we
obtain in this study, where we measure $\alpha<0$ for the various CO tracers
in the n30 and n100 model, given our fitting range. However, this is primarily
due to the fact that our $\Delta$-variance spectra peak at roughly the scale
of small carbon monoxide clumps in the low-density clouds (see Section 4.2).
Accordingly, we find more cloud structures on the scales of those localized CO
structures, leading to negative $\alpha$ slopes in Figure 5. A similar effect
can be seen in the $\Delta$-variance analysis of 13CO ($J=1\rightarrow 0$)
maps of the outer Galaxy shown in Figure 3 of Stutzki et al. (1998). On small
scales, a positive power-law slope is recovered, but on large scales $\alpha$
becomes negative as one starts to smooth on scales larger than the typical
size of the observed MCs. The same situation is also illustrated in the
$\Delta$-variance spectra in Schneider et al. (2011) for the 13CO
($J=1\rightarrow 0$) molecular line survey of Cygnus X as well as in a 13CO
($J=2\rightarrow 1$) study of the Perseus cloud in Sun et al. (2006), which
peak at intermediate cloud scales. However, a $\Delta$-variance analysis of
the small sub-regions NGC1333 and L1455 in the Perseus molecular cloud complex
shows that $\alpha$ is always positive (see Fig. 6 in Sun et al., 2006), which
is in agreement with our finding that a positive slope $\alpha$ might be
related to a $\Delta$-variance analysis at lower spatial cloud scales, probing
only very localized CO structures. Moreover, Alves de Oliveira et al. (2014)
analysed the Chamaeleon molecular cloud complex using data from the Herschel
Gould Belt Survey. This complex encompasses three MCs with different star
formation histories. One of them is a quiescent cloud, which should be best
comparable to our numerical simulations, since we neglect the effects of self-
gravity. The $\Delta$-variance analysis of this MC also shows a clear break at
mid to large scales (see cloud Cha III in Fig. 5 in Alves de Oliveira et al.,
2014), also displaying a characteristic spatial scale in the column density
structure.

In our high-density n300 run, the CO gas is distributed over the whole MC and
not only confined to small dense fragments. In this case, we also find cloud
structures on larger spatial scales, leading to positive $\alpha$ slopes in
Figure 5.

###  4.6 Model limitations

Since we are running numerical simulations, we have to keep in mind that our
runs are subject to various physical simplifications. In the current analysis,
we focus on the impact of the chemistry and the opacity of the gas on our
results, thus neglecting other physical processes that could bias our
analysis. In first instance, we do not account for self-gravity and thus also
do not model star formation or stellar feedback. More specifically we do not
account for stellar radiation, SN feedback and other physical processes. We
also do not include any large-scale dynamics, e.g. spiral arms or galactic
rotation, although these are unlikely to be important on the 20 pc scale
studied here. Nevertheless, we can infer useful information about how the
chemical composition of the gas and the opacity affect the $\Delta$-variance
analysis. We also note that our results depend only weakly on the resolution
and that this concerns mostly CO as a tracer molecule, as we show in Appendix
A. For future investigations, we want to analyze simulations that span a wider
range of physical parameters, e.g. with different levels of the external
radiation field, varying metallicities or additional physical processes, in
order to find out how they affect the statistics.

##  5 Summary and Conclusions

We analyzed $\Delta$-variance spectra of MCs with time-dependent chemistry and
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 12CO ($J=1\rightarrow 0$) and
13CO ($J=1\rightarrow 0$) intensity (both with radiative transfer). In each
case, we computed $\Delta$-variance spectra for maps of centroid velocities
(CV), integrated intensities and column densities and analyzed the structural
behavior of MCs in numerical simulations. We report the following findings:

  * •

We compute $\Delta$-variance spectra of maps of centroid velocities and fit a
power-law $\sigma_{\Delta}^{2}(\ell)\propto\ell^{\alpha}$, in order to
characterize the properties of the turbulent hierarchy in the MCs. This power-
law can be translated into a linewidth-size relation, i.e.
$\sigma_{\Delta}(\ell)\propto\ell^{\gamma}$ with $\gamma=\alpha/2$, readily
comparable to slopes that can be derived from spectral observations. We find
the slopes $\alpha$ of both the total and H2 density models to be
significantly steeper than the slopes of the different CO tracers, which
underestimate the former by a factor of $\sim 1.5-3.0$ (see Section 3.1).

  * •

The slopes $\alpha$ derived from the CV maps range from 0.8 to 1.3 for the
total and H2 density, while $\alpha$ for the various CO tracers range from 0.3
to 0.8 (see Section 3.1). However, we also note that the specific choice of
the fitting range might cause further variations of the slopes by $\pm 0.1$.

  * •

Although we find slight variations between the different slopes $\alpha$ for
our various CO models obtained from the CV maps, the impact of the optical
depth effects on the spectra computed on maps of centroid velocities remains
limited (see Table 2).

  * •

This is different in the case of the integrated intensity and column density.
The $\Delta$-variance computed from these maps is strongly affected by optical
depth effects. The CO tracers exhibit a very different spatial scaling
behavior compared to the total and H2 density models (see Table 3 and Figure
5).

  * •

We report a critical number density threshold of $\sim 100\,$cm-3 at which the
spectral slopes $\alpha$ of the CO tracers switch sign for the
$\Delta$-variance of integrated intensity and column density maps. We conclude
that carbon monoxide traces the total cloud structure well only if the average
cloud density lies significantly above this threshold. If the mean density in
the cloud is significantly smaller than this limit, the observable CO gas does
not properly trace the statistical properties of the H2 gas in the cloud (see
Section 3.2).

  * •

The $\Delta$-variance spectra computed on maps of integrated intensity and
column density provide a useful statistical measure in order to infer
important information about the distribution of gas within a cloud. We also
argue that peaks in the $\Delta$-variance spectra correspond to characteristic
scales of the morphological structure of the system (see Section 3.2).

  * •

Our findings are consistent with previous $\Delta$-variance studies using CO
line observations or measurements of the continuum (see Section 4.4 and 4.5).

## Acknowledgements

We thank Volker Ossenkopf for informative discussions about the usage of the
$\Delta$-variance and for providing his IDL routines in order to compute the
$\Delta$-variance spectra. We also thank Lukas Konstandin for stimulating
discussions about the theory of turbulence as well as the referee for a very
constructive and detailed report, which helped to improve the paper. EB, SCOG
and RSK acknowledge support from the Deutsche Forschungsgemeinschaft (DFG) via
the SFB 881 (sub-projects B1, B2, B5 and B8) “The Milky Way System”, and the
SPP (priority program) 1573, “Physics of the ISM”. Furthermore, EB
acknowledges financial support from the Konrad-Adenauer-Stiftung (KAS) via
their “Promotionsförderung”. Some of the simulations presented 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 Research Council under the European
Community’s Seventh Framework Programme (FP7/2007-2013) via the ERC Advanced
Grant ”STARLIGHT: Formation of the First Stars” (project number 339177).

##  Appendix A Resolution study

Figure 6: Same as Fig. 4, but with runs of different resolutions with
$512^{3}$ and $256^{3}$ grid cells with a fixed initial number density of
$n_{0}=100\,$cm-3, computed for the centroid velocity (CV) maps. We find
similar $\alpha$ values for the total density and the H2 density models, while
the values $\alpha$ differ by up to $\sim 20-30\%$ between the different
resolution models for the various CO tracers.

We study the influence of the numerical resolution on the results of the
$\Delta$-variance. Therefore, we have performed runs with resolutions of
$512^{3}$ and $256^{3}$ grid cells and evaluate the spectra and slopes for the
CV maps for all chemical components for a fixed initial number density of
$n_{0}=100\,$cm-3. The results and their interpretation is the same for all
other density models. Fig. 6 shows spectra of the $\Delta$-variance with the
corresponding slopes $\alpha$, in analogy to Fig. 4. The fitting range for the
$256^{3}$ model is downscaled by a factor of 2 compared to the $512^{3}$
model, i.e. we fit from 25 to 64 cells in the spatial domain. We find similar
$\alpha$ values for the total density and the H2 density models. However, the
slope values $\alpha$ differ by up to $\sim 20-30\%$ between the different
resolution models for the various CO tracers. This is because CO is mainly
located in dense regions of the cloud (Bertram et al., 2014), which can be
resolved more accurately at a higher resolution, leading to significant
differences between the two resolution models. These results agree with the
results in the resolution study of the structure function analysis presented
in Bertram et al. (2015a). However, the variations caused by the effect of
resolution are rather small, since $\alpha$ also strongly depends on the
specific choice of the fitting range, which might also cause slope variations
by about $\pm 0.1$. Nevertheless, we find a similar relative scaling behavior
in the two resolution models between the spectra of the different chemical
components.

##  Appendix B Comparison of spectra with different filter functions

Figure 7: Same as Fig. 4, but with the $\Delta$-variance spectra for the H2
and CO density models for a fixed initial number density of $n_{0}=100\,$cm-3
and a resolution of $512^{3}$ grid cells, computed for the centroid velocity
(CV) maps. In order to analyze the impact of the filter function and its
diameter ratio on our results, we compute the H2 and CO density spectra using
both a Mexican hat with a diameter ratio of 1.5 and a French hat with a
diameter ratio of 3.0. We do not find any significant differences in the
slopes $\alpha$ within the fitting errors if we use a another filter function
for the $\Delta$-variance analysis combined with a different diameter ratio.

We also study the influence of the filter function and the choice of the
specific diameter ratio on our $\Delta$-variance spectra. Therefore, as an
example, Fig. 7 shows the $\Delta$-variance spectra for the H2 and CO density
models for a fixed initial number density of $n_{0}=100\,$cm-3 and a
resolution of $512^{3}$ grid cells, computed on maps of centroid velocities.
For each model, we evaluate the spectra using different filter functions and
diameter ratios. In particular, we compute the $\Delta$-variance spectra using
a Mexican hat with a diameter ratio of 1.5 as well as a French hat with a
diameter ratio of 3.0. In analogy to Section 3.1, we fit a power-law within a
given fitting range to the spectra and compare the slopes with each other,
which are shown in Fig. 7. Thereby, we do not find any significant differences
between the slopes derived from spectra with various filter functions and
diameter ratios within the fitting errors for our models. The individual
normalizations of the spectra are caused by the variable choice of the
diameter ratio, affecting the computation of the $\Delta$-variance analysis
(see, e.g. Ossenkopf et al., 2008a). However, the shape of the individual
spectra for one chemical model over various spatial scales is also
approximately conserved.

## References

  * Alves de Oliveira et al. (2014) Alves de Oliveira C., Schneider N., Merín B., Prusti T., Ribas Á., Cox N. L. J., Vavrek R., Könyves V., Arzoumanian D., Puga E., Pilbratt G. L., Kóspál Á., André P. e. a., 2014, A&A, 568, A98
  * André et al. (2010) André P., Men’shchikov A., Bontemps S., Könyves V., Motte F., Schneider N., Didelon P., Minier V., Saraceno P., Ward-Thompson D., di Francesco J., White G., Molinari S., Testi L., Abergel 2010, A&A, 518, L102
  * Bensch et al. (2001) Bensch F., Stutzki J., Ossenkopf V., 2001, A&A, 366, 636
  * Benzi et al. (1993) Benzi R., Ciliberto S., Tripiccione R., Baudet C., Massaioli F., Succi S., 1993, Phys. Rev. E, 48, R29
  * Bertram et al. (015b) Bertram E., Glover S. C. O., Clark P., Klessen R. S., 2015b, MNRAS
  * Bertram et al. (015a) Bertram E., Konstandin L., Shetty R., Glover S. C. O., Klessen R. S., 2015a, MNRAS, 446, 3777
  * Bertram et al. (2014) Bertram E., Shetty R., Glover S. C. O., Klessen R. S., Roman-Duval J., Federrath C., 2014, MNRAS, 440, 465
  * Boldyrev et al. (2002) Boldyrev S., Nordlund Å., Padoan P., 2002, Physical Review Letters, 89, 031102
  * Brunt & Heyer (2002a) Brunt C. M., Heyer M. H., 2002a, ApJ, 566, 276
  * Brunt & Heyer (2002b) Brunt C. M., Heyer M. H., 2002b, ApJ, 566, 289
  * Burgers (1948) Burgers J. M., 1948, Adv. Appl. Mech., 1, 171
  * Burkhart et al. (2014) Burkhart B., Lazarian A., Leão I. C., de Medeiros J. R., Esquivel A., 2014, ApJ, 790, 130
  * Burkhart et al. (013a) Burkhart B., Lazarian A., Ossenkopf V., Stutzki J., 2013a, ApJ, 771, 123
  * Burkhart et al. (013b) Burkhart B., Ossenkopf V., Lazarian A., Stutzki J., 2013b, ApJ, 771, 122
  * Draine (1978) Draine B. T., 1978, ApJS, 36, 595
  * Dubrulle (1994) Dubrulle B., 1994, Phys. Rev. Lett., 73, 959
  * Dullemond (2012) Dullemond C. P., , 2012, RADMC-3D: A multi-purpose radiative transfer tool, Astrophysics Source Code Library
  * Elia et al. (2014) Elia D., Strafella F., Schneider N., Paladini R., Vavrek R., Maruccia Y., Molinari S., Noriega-Crespo A., Pezzuto S., Rygl K. L. J., Di Giorgio A. M., Traficante A., Schisano E., Calzoletti L. e. a., 2014, ApJ, 788, 3
  * Elmegreen & Scalo (2004) Elmegreen B. G., Scalo J., 2004, ARA&A, 42, 211
  * Esquivel & Lazarian (2005) Esquivel A., Lazarian A., 2005, ApJ, 631, 320
  * Esquivel et al. (2007) Esquivel A., Lazarian A., Horibe S., Cho J., Ossenkopf V., Stutzki J., 2007, MNRAS, 381, 1733
  * Falgarone et al. (1998) Falgarone E., Panis J.-F., Heithausen A., Perault M., Stutzki J., Puget J.-L., Bensch F., 1998, A&A, 331, 669
  * Federrath et al. (2009) Federrath C., Klessen R. S., Schmidt W., 2009, ApJ, 692, 364
  * Federrath et al. (2010) Federrath C., Roman-Duval J., Klessen R. S., Schmidt W., Mac Low M.-M., 2010, A&A, 512, A81
  * Frisch & Kolmogorov (1995) Frisch U., Kolmogorov A., 1995, Turbulence: The Legacy of A.N. Kolmogorov. Cambridge University Press
  * Glover & Clark (2012) Glover S. C. O., Clark P. C., 2012, MNRAS, 421, 116
  * Glover et al. (2010) Glover S. C. O., Federrath C., Mac Low M.-M., Klessen R. S., 2010, MNRAS, 404, 2
  * Glover & Mac Low (2007) Glover S. C. O., Mac Low M.-M., 2007, ApJS, 169, 239
  * Glover & Mac Low (2011) Glover S. C. O., Mac Low M.-M., 2011, MNRAS, 412, 337
  * Habing (1968) Habing H. J., 1968, Bull. Astron. Inst. Netherlands, 19, 421
  * Hayes et al. (2006) Hayes J. C., Norman M. L., Fiedler R. A., Bordner J. O., Li P. S., Clark S. E., ud-Doula A., Mac Low M.-M., 2006, ApJS, 165, 188
  * Heithausen & Thaddeus (1990) Heithausen A., Thaddeus P., 1990, ApJ, 353, L49
  * Hennebelle & Falgarone (2012) Hennebelle P., Falgarone E., 2012, A&A Rev., 20, 55
  * Heyer & Brunt (2004) Heyer M. H., Brunt C. M., 2004, ApJ, 615, L45
  * Heyer & Schloerb (1997) Heyer M. H., Schloerb F. P., 1997, ApJ, 475, 173
  * Heyer et al. (2006) Heyer M. H., Williams J. P., Brunt C. M., 2006, ApJ, 643, 956
  * Hily-Blant et al. (2008) Hily-Blant P., Falgarone E., Pety J., 2008, A&A, 481, 367
  * Johnston et al. (2014) Johnston K. G., Beuther H., Linz H., Schmiedeke A., Ragan S. E., Henning T., 2014, A&A, 568, A56
  * Klessen (2000) Klessen R. S., 2000, ApJ, 535, 869
  * Klessen & Glover (2014) Klessen R. S., Glover S. C. O., 2014, ArXiv e-prints, 1412.5182
  * Klessen et al. (2000) Klessen R. S., Heitsch F., Mac Low M.-M., 2000, ApJ, 535, 887
  * Kolmogorov (1941) Kolmogorov A., 1941, Akademiia Nauk SSSR Doklady, 30, 301
  * Konstandin et al. (2012) Konstandin L., Federrath C., Klessen R. S., Schmidt W., 2012, Journal of Fluid Mechanics, 692, 183
  * Konstandin et al. (2015) Konstandin L., Shetty R., Girichidis P., Klessen R. S., 2015, MNRAS, 446, 1775
  * Kritsuk et al. (2007) Kritsuk A. G., Norman M. L., Padoan P., Wagner R., 2007, ApJ, 665, 416
  * Larson (1981) Larson R. B., 1981, MNRAS, 194, 809
  * Lazarian & Esquivel (2003) Lazarian A., Esquivel A., 2003, ApJ, 592, L37
  * Lazarian & Pogosyan (2000) Lazarian A., Pogosyan D., 2000, ApJ, 537, 720
  * Lazarian & Pogosyan (2004) Lazarian A., Pogosyan D., 2004, ApJ, 616, 943
  * Lis et al. (1996) Lis D. C., Pety J., Phillips T. G., Falgarone E., 1996, ApJ, 463, 623
  * Mac Low (1999) Mac Low M.-M., 1999, ApJ, 524, 169
  * Mac Low & Klessen (2004) Mac Low M.-M., Klessen R. S., 2004, Reviews of Modern Physics, 76, 125
  * Mac Low et al. (1998) Mac Low M.-M., Klessen R. S., Burkert A., Smith M. D., 1998, Physical Review Letters, 80, 2754
  * Mac Low & Ossenkopf (2000) Mac Low M.-M., Ossenkopf V., 2000, A&A, 353, 339
  * McKee & Ostriker (2007) McKee C. F., Ostriker E. C., 2007, ARA&A, 45, 565
  * Meyerdierks & Heithausen (1996) Meyerdierks H., Heithausen A., 1996, A&A, 313, 929
  * Miville-Deschênes et al. (2010) Miville-Deschênes M.-A., Martin P. G., Abergel A., Bernard J.-P., Boulanger F., Lagache G., Anderson L. D., André P., Arab H., Baluteau J.-P., Blagrave K., Bontemps S., Cohen 2010, A&A, 518, L104
  * Nelson & Langer (1999) Nelson R. P., Langer W. D., 1999, ApJ, 524, 923
  * Norman (2000) Norman M. L., 2000, in Arthur S. J., Brickhouse N. S., Franco J., eds, Revista Mexicana de Astronomia y Astrofisica Conference Series Vol. 9 of Revista Mexicana de Astronomia y Astrofisica, vol. 27, Introducing ZEUS-MP: A 3D, Parallel, Multiphysics Code for Astrophysical Fluid Dynamics. pp 66–71
  * Ossenkopf et al. (1998) Ossenkopf V., Bensch F., Zielinsky M., 1998, Interstellar Turbulence. Cambridge, 252
  * Ossenkopf et al. (2006) Ossenkopf V., Esquivel A., Lazarian A., Stutzki J., 2006, A&A, 452, 223
  * Ossenkopf et al. (2001) Ossenkopf V., Klessen R. S., Heitsch F., 2001, A&A, 379, 1005
  * Ossenkopf et al. (2008a) Ossenkopf V., Krips M., Stutzki J., 2008a, A&A, 485, 917
  * Ossenkopf et al. (2008b) Ossenkopf V., Krips M., Stutzki J., 2008b, A&A, 485, 719
  * Ossenkopf & Mac Low (2002) Ossenkopf V., Mac Low M.-M., 2002, A&A, 390, 307
  * Röllig & Ossenkopf (2013) Röllig M., Ossenkopf V., 2013, A&A, 550, A56
  * Roman-Duval et al. (2011) Roman-Duval J., Federrath C., Brunt C., Heyer M., Jackson J., Klessen R. S., 2011, ApJ, 740, 120
  * Roman-Duval et al. (2010) Roman-Duval J., Jackson J. M., Heyer M., Rathborne J., Simon R., 2010, ApJ, 723, 492
  * Rosolowsky et al. (1999) Rosolowsky E. W., Goodman A. A., Wilner D. J., Williams J. P., 1999, ApJ, 524, 887
  * Russeil et al. (2013) Russeil D., Schneider N., Anderson L. D., Zavagno A., Molinari S., Persi P., Bontemps S., Motte F., Ossenkopf V., André P., Arzoumanian D., Bernard J.-P., Deharveng L., Didelon P. e. a., 2013, A&A, 554, A42
  * Scalo & Elmegreen (2004) Scalo J., Elmegreen B. G., 2004, ARA&A, 42, 275
  * Schneider et al. (2011) Schneider N., Bontemps S., Simon R., Ossenkopf V., Federrath C., Klessen R. S., Motte F., André P., Stutzki J., Brunt C., 2011, A&A, 529, A1
  * Sembach et al. (2000) Sembach K. R., Howk J. C., Ryans R. S. I., Keenan F. P., 2000, ApJ, 528, 310
  * She & Leveque (1994) She Z.-S., Leveque E., 1994, Phys. Rev. Lett., 72, 336
  * Shetty et al. (011a) Shetty R., Glover S. C., Dullemond C. P., Klessen R. S., 2011a, MNRAS, 412, 1686
  * Shetty et al. (011b) Shetty R., Glover S. C., Dullemond C. P., Ostriker E. C., Harris A. I., Klessen R. S., 2011b, MNRAS, 415, 3253
  * Sobolev (1957) Sobolev V. V., 1957, SvA, 1, 678
  * Solomon et al. (1987) Solomon P. M., Rivolo A. R., Barrett J., Yahil A., 1987, ApJ, 319, 730
  * Stutzki et al. (1998) Stutzki J., Bensch F., Heithausen A., Ossenkopf V., Zielinsky M., 1998, A&A, 336, 697
  * Sun et al. (2006) Sun K., Kramer C., Ossenkopf V., Bensch F., Stutzki J., Miller M., 2006, A&A, 451, 539
  * Szűcs et al. (2014) Szűcs L., Glover S. C. O., Klessen R. S., 2014, MNRAS, 445, 4055
  * Visser et al. (2009) Visser R., van Dishoeck E. F., Black J. H., 2009, A&A, 503, 323
  * Watson et al. (1976) Watson W. D., Anicich V. G., Huntress Jr. W. T., 1976, ApJ, 205, L165