Comparing the statistics of interstellar turbulence in simulations and observations
Solenoidal versus compressive turbulence forcing^{}
C. Federrath^{1,2,3}  J. RomanDuval^{4,5}  R. S. Klessen^{1}  W. Schmidt^{6,7}  M.M. Mac Low^{2,3}
1  Zentrum für Astronomie der Universität Heidelberg, Institut für
Theoretische Astrophysik, AlbertUeberleStr. 2, 69120
Heidelberg, Germany
2  MaxPlanckInstitut für Astronomie, Königstuhl 17, 69117
Heidelberg, Germany
3  Department of Astrophysics, American Museum of Natural History,
Central Park West at 79th Street, New York, NY 100245192, USA
4  Space Telescope Science Institute, 3700 San Martin Drive,
Baltimore, MD 21218, USA
5  Astronomy Department, Boston University, 725 Commonwealth Avenue,
Boston, MA 02215, USA
6  Institut für Astrophysik, Universität Göttingen,
FriedrichHundPlatz 1, 37077 Göttingen, Germany
7  Lehrstuhl für Astronomie, Institut für Theoretische Physik und
Astrophysik, Universität Würzburg, Am Hubland, 97074 Würzburg, Germany
Received 7 May 2009 / Accepted 15 January 2010
Abstract
Context. Density and velocity fluctuations on
virtually all scales observed with modern telescopes show that
molecular clouds (MCs) are turbulent. The forcing and structural
characteristics of this turbulence are, however, still poorly
understood.
Aims. To shed light on this subject, we study two
limiting cases of turbulence forcing in numerical experiments:
solenoidal (divergencefree) forcing and compressive (curlfree)
forcing, and compare our results to observations.
Methods. We solve the equations of hydrodynamics on
grids with up to 1024^{3} cells for purely
solenoidal and purely compressive forcing. Eleven lowerresolution
models with different forcing mixtures are also analysed.
Results. Using Fourier spectra and variance,
we find velocity dispersionsize relations consistent with observations
and independent numerical simulations, irrespective of the type of
forcing. However, compressive forcing yields stronger compression at
the same rms Mach number than solenoidal forcing, resulting in a three
times larger standard deviation of volumetric and column density
probability distributions (PDFs). We compare our results to different
characterisations of several observed regions, and find evidence of
different forcing functions. Column density PDFs in the
Perseus MC suggest the presence of a mainly
compressive forcing agent within a shell, driven by a massive star.
Although the PDFs are close to lognormal, they have nonGaussian
skewness and kurtosis caused by intermittency. Centroid velocity
increments measured in the Polaris Flare on intermediate
scales agree with solenoidal forcing on that scale. However, variance
analysis of the column density in the Polaris Flare suggests
that turbulence is driven on large scales, with a significant
compressive component on the forcing scale. This indicates that,
although likely driven with mostly compressive modes on large scales,
turbulence can behave like solenoidal turbulence on smaller scales.
Principal component analysis of G2162.5 and most of the
Rosette MC agree with solenoidal forcing, but the
interior of an ionised shell within the Rosette MC
displays clear signatures of compressive forcing.
Conclusions. The strong dependence of the density
PDF on the type of forcing must be taken into account in any theory
using the PDF to predict properties of star formation. We supply a
quantitative description of this dependence. We find that different
observed regions show evidence of different mixtures of compressive and
solenoidal forcing, with more compressive forcing occurring primarily
in sweptup shells. Finally, we emphasise the role of the sonic scale
for protostellar core formation, because core formation close to the
sonic scale would naturally explain the observed subsonic velocity
dispersions of protostellar cores.
Key words: hydrodynamics  ISM: clouds  ISM: kinematics and dynamics  methods: numerical  methods: statistical  turbulence
1 Introduction
Studying the density and velocity distributions of interstellar gas provides essential information about virtually all physical processes relevant to the dynamical evolution of the interstellar medium (ISM). Along with gravity, magnetic fields and the thermodynamics of the gas, supersonic turbulence plays a fundamental role in determining the density and velocity statistics of the ISM (e.g., Scalo et al. 1998). Thus, supersonic turbulence is considered a key process for star formation (Scalo & Elmegreen 2004; Elmegreen & Scalo 2004; McKee & Ostriker 2007; Mac Low & Klessen 2004).
In this paper, we continue our analysis of the density probability distribution function (PDF) obtained in numerical experiments of driven supersonic isothermal turbulence. Understanding the density PDF and its turbulent origin is essential, because it is a key ingredient for analytical models of star formation: the turbulent density PDF is used to explain the stellar initial mass function (Hennebelle & Chabrier 2009,2008; Padoan & Nordlund 2002), the star formation rate (Padoan & Nordlund 2009; Krumholz & McKee 2005; Krumholz et al. 2009), the star formation efficiency (Elmegreen 2008), and the KennicuttSchmidt relation on galactic scales (Elmegreen 2002; Tassis 2007; Kravtsov 2003). In Federrath et al. (2008b), we found that supersonic turbulence driven by a purely compressive (curlfree) force field yields a density PDF with roughly three times larger standard deviation compared to solenoidal (divergencefree) turbulence forcing, which strongly affects the results obtained in these analytical models. Here, we want to compare our results for the density PDF to observations of column density PDFs (e.g., Goodman et al. 2009).
Moreover, in Federrath et al. (2009) we investigated the fractal density distribution of our two models with solenoidal and compressive turbulence forcing, which showed that compressive forcing yields a significantly lower fractal dimension ( ) compared to solenoidal forcing ( ). In the present contribution, we consider the scaling of centroid velocity increments computed for these models, and we compare them to observations of the Polaris Flare by HilyBlant et al. (2008). We additionally used principal component analysis and compared our results to observations of the G2162.5 (Maddalena's Cloud) and the Rosette MC by Heyer et al. (2006).
Our results indicate that interstellar turbulence is driven by mixtures of solenoidal and compressive forcing. The ratio between solenoidal and compressive modes of the turbulence forcing may vary strongly across different regions of the ISM. This provides an explanation for the apparent lack of correlation between turbulent density and velocity dispersions found in observations (e.g., Pineda et al. 2008; Goodman et al. 2009). We conclude that solenoidal forcing is more likely to be realised in quiescent regions with low star formation activity as in the Polaris Flare and in Maddalena's Cloud. On the other hand, in regions of enhanced stellar feedback, compressive forcing leads to larger standard deviations of the density PDFs, as seen in one of the subregions of the Perseus MC surrounding a central B star. Moreover, compressive forcing exhibits a higher scaling exponent of principal component analysis than solenoidal forcing. This higher scaling exponent is consistent with the measured scaling exponent for the interior of an ionising shell in the Rosette MC.
In Sect. 2, we explain the numerical setup and turbulence forcing used for the present study. We discuss our results obtained using PDFs, centroid velocity increments, principal component analysis, Fourier spectrum functions, and variance analyses in Sects. 37, respectively. In each of these sections, we compare the turbulence statistics obtained for solenoidal and compressive forcing with observational data available in the literature. In Sect. 8, we discuss the possibility that transonic prestellar cores typically form close to the sonic scale in a globally supersonic, turbulent medium. Section 9 provides a list of the limitations in our comparison of numerical simulations with observations. A summary of our results and conclusions is given in Sect. 10.
2 Simulations and methods
The piecewise parabolic method (Colella
& Woodward 1984), implemented in the astrophysical
code FLASH3 (Dubey
et al. 2008; Fryxell et al. 2000)
was used to integrate the equations of hydrodynamics on
threedimensional (3D) periodic uniform grids with 256^{3},
512^{3}, and 1024^{3} grid points.
Since isothermal gas is assumed throughout this study, it is
convenient to define
as the natural logarithm of the density divided by the mean density in the system. For isothermal gas, the pressure, , is proportional to the density with the constant sound speed . The equations of hydrodynamics solved here are consequently given by
where denotes the velocity of the gas. An energy equation is not needed, because the gas is isothermal. The assumption of isothermal gas is very crude, but may still provide an adequate physical approximation to the real thermodynamics in dense molecular gas (Wolfire et al. 1995; Pavlovski et al. 2006). We discuss further limitations of our simulations in Sect. 9. The stochastic forcing term is used to drive turbulent motions.
2.1 Forcing module
Equations (2) and (3) have been solved before in the context of molecular cloud dynamics, studying compressible turbulence with either solenoidal (divergencefree) forcing or with a 2:1 mixture of solenoidal to compressive modes in the turbulence forcing (e.g., Dib et al. 2008; Stone et al. 1998; Padoan et al. 2004; Mac Low et al. 1998; Offner et al. 2008; Klessen et al. 2000; Mac Low 1999; Padoan et al. 1997; Klessen 2001; Li et al. 2003; Schmidt et al. 2009; BallesterosParedes et al. 2006; Jappsen et al. 2005; Boldyrev et al. 2002; Heitsch et al. 2001; Kritsuk et al. 2007; Kissmann et al. 2008). The case of a 2:1 mixture of solenoidal to compressive modes is the natural result obtained for 3D forcing, if no Helmholtz decomposition (see below) is performed. Then, the solenoidal modes occupy two of the three available spatial dimensions on average, while the compressive modes only occupy one (Federrath et al. 2008b; Elmegreen & Scalo 2004). In the present study, the solenoidal forcing case is thus also used as a control run for comparison with previous studies using solenoidal forcing. However, we additionally applied purely compressive (curlfree) forcing and analysed the resulting turbulence statistics in detail. Each simulation at a resolution of 1024^{3} grid cells consumed roughly . Therefore, we concentrated on two extreme cases of turbulence forcing with high resolution: (1) the widely adopted purely solenoidal forcing ( ), and (2) purely compressive forcing ( ). However, we also studied eleven simulations at numerical resolution of 256^{3} in which we smoothly varied the forcing from purely solenoidal to purely compressive by producing eleven different forcing mixtures.
The forcing term is often modelled with a spatially static pattern, for which the amplitude is adjusted in time following the methods introduced by Mac Low et al. (1998) and Stone et al. (1998). This results in a roughly constant energy input on large scales. Other studies model the random forcing term such that it can vary in time and space (e.g., Federrath et al. 2008b; Schmidt et al. 2009; Padoan et al. 2004; Kritsuk et al. 2007). Here, we used the OrnsteinUhlenbeck (OU) process to model , which belongs to the latter type. The OU process is a welldefined stochastic process with a finite autocorrelation timescale. It can be used to excite turbulent motions in 3D, 2D, and 1D simulations as explained in Eswaran & Pope (1988) and Schmidt et al. (2006). Using an OU process enables us to control the autocorrelation timescale T of the forcing. The concept of using the OU process to excite turbulence and the projections in Fourier space necessary to get solenoidal and compressive force fields are described below.
The OU process is a stochastic differential equation
describing the evolution of the forcing term
in Fourier space (kspace):
The first term on the right hand side is a diffusion term. This term is modelled using a Wiener process , which adds a Gaussian random increment to the vector field given in the previous time step . Wiener processes are random processes, such that
(5) 
where denotes the 3D, 2D, or 1D version of a Gaussian distribution with zero mean and standard deviation . This is followed by a projection with the projection tensor in Fourier space. In index notation, the projection operator reads
where is the Kronecker symbol, and and are the fully solenoidal and the fully compressive projection operators, respectively. The projection operator serves to construct a purely solenoidal force field by setting . For , a purely compressive force field is obtained. Any combination of solenoidal and compressive modes can be constructed by choosing . By changing the parameter , we can thus set the power of compressive modes with respect to the total power of the forcing. The analytical ratio of compressive power to total power can be derived from Eq. (6) by evaluating the norm of the compressive component of the projection tensor,
and by evaluating the norm of the full projection tensor
The result of the last equation depends on the dimensionality D=1,2,3 of the forcing, because the norm of the Kronecker symbol and 3 in one, two and three dimensions, respectively. The ratio of Eqs. (7) and (8) gives the ratio of compressive forcing power to the total forcing power as a function of the parameter :
Figure 1 provides a graphical representation of this ratio for the 1D, 2D, and 3D case. For comparison, we plot numerical values of the forcing ratio obtained in eleven 3D and 2D hydrodynamical runs with resolutions of 256^{3} and 1024^{2} grid points, in which we have varied the forcing parameter from purely compressive forcing () to purely solenoidal forcing () in the range , separated by . Note that a natural mixture of forcing modes is obtained for , which leads to for 3D turbulence, and for 2D turbulence. A simple way to understand this natural ratio is to consider longitudinal and transverse waves. In 3D, the longitudinal waves occupy one of the three spatial dimensions, while the transverse waves occupy two of the three on average. Thus, the longitudinal (compressive) part has a power of 1/3, while the transverse (solenoidal) part has a power of 2/3 in 3D. In 2D, the natural ratio is 1/2, because longitudinal and transverse waves are evenly distributed in two dimensions.
Figure 1: Ratio of compressive power to the total power in the turbulence force field. The solid lines labelled with 1D, 2D, and 3D show the analytical expectation for this ratio, Eq. (9), as a function of the forcing parameter for one, two and threedimensional forcing, respectively. The diamonds and squares show results of numerical simulations in 3D and 2D with , separated by . Those models were run at a numerical resolution of 256^{3} and 1024^{2} grid points in 3D and 2D, respectively. The two extreme forcing cases of purely solenoidal forcing () and purely compressive forcing () are indicated as ``sol'' and ``comp'', respectively. Note that in any 1D model, all power is in the compressive component, and thus , independent of . 

Open with DEXTER 
The second term on the right hand side of Eq. (4) is a drift term, which models the exponentially decaying correlation of the force field with itself. Thus, the autocorrelation timescale of the forcing is denoted by T. We set the autocorrelation timescale equal to the dynamical timescale T=L/(2V) on the scale of energy injection, where L is the size of the computational domain, and is the rms Mach number in all runs. The autocorrelation timescale is therefore equal to the decay time constant in supersonic hydrodynamic and magnetohydrodynamic turbulence driven on large scales (Stone et al. 1998; Mac Low 1999). The forcing amplitude is a paraboloid in 3D Fourier space, only containing power on the largest scales in a small interval of wavenumbers peaking at k=2, which corresponds to half of the box size L/2. The effects of varying the scale of energy input were investigated by Mac Low (1999), Klessen et al. (2000), Heitsch et al. (2001) and VázquezSemadeni et al. (2003). Here, we consider largescale stochastic forcing, which is closer to the observational data (e.g., Ossenkopf & Mac Low 2002; Brunt et al. 2009). This type of forcing models the kinetic energy input from largescale turbulent fluctuations, breaking up into smaller structures. Kinetic energy cascades down to smaller and smaller scales, and thereby effectively drives turbulent fluctuations on scales smaller than the turbulence injection scale.
We have verified that our results are not sensitive to the general approach of using an OrnsteinUhlenbeck process for the turbulence forcing. For instance, we have used an almost static forcing pattern, which is obtained in the limit in test simulations. We have furthermore checked that the particular choice of Fourier amplitudes did not affect our results by using a band spectrum instead of a parabolic forcing spectrum. Varying these parameters did not strongly affect our results. In contrast, changing from (solenoidal forcing) to (compressive forcing) always led to significant changes in the turbulence statistics.
2.2 Initial conditions and postprocessing
Starting from a uniform density distribution and zero velocities, the forcing excites turbulent motions. Equations (2) and (3) have been evolved for ten dynamical times T, which allows us to study a large sample of realisations of the turbulent flow. Compressible turbulence reached a statistically invariant state within 2 T (Federrath et al. 2009). This allows us to average all statistical measures over 8 T separated by 0.1 T in the fully developed regime. We are thus able to average over 81 different realisations of the turbulence to improve statistical significance. The 1 temporal fluctuations obtained from this averaging procedure are indicated as error bars for the PDFs, centroid velocity increments, principal component analysis, Fourier spectra and variance analyses in the following sections and in all figures showing error bars throughout this study. The forcing amplitude was adjusted to excite a turbulent flow with an rms Mach number in all cases. We use the rms Mach number as the control parameter, because this dimensionless number determines most of the physical properties of scaleinvariant turbulent flows and is often used to derive important flow statistics such as the standard deviation of the density distribution. However, in the next section we show that the latter depends sensitively on the turbulence forcing parameter as well.
Figure 2 (top panels) shows column density fields projected along the zaxis from a randomly selected snapshot at time t=2 T in the regime of fully developed, statistically stationary turbulence for solenoidal (left) versus compressive forcing (right). This regime was reached after 2 dynamical times T, which is shown in Fig. 3 for the minimum and maximum logarithmic densities s (top panel) and rms curl and divergence of the velocity field (bottom panel) as a function of the dynamical time. It is evident that compressive forcing produces higher density contrasts, resulting in higher density peaks and bigger voids compared to solenoidal forcing.
Figure 2: Maps showing density ( top), vorticity ( middle) and divergence ( bottom) in projection along the zaxis at time t=2 T as an example for the regime of statistically fully developed, compressible turbulence for solenoidal forcing ( left) and compressive forcing ( right). Top panels: column density fields in units of the mean column density. Both maps show three orders of magnitude in column density with the same scaling and magnitudes for direct comparison. Middle panels: projections of the modulus of the vorticity . Regions of intense vorticity appear to be elongated filamentary structures often coinciding with positions of intersecting shocks. Bottom panels: projections of the divergence of the velocity field showing the positions of shocks. Negative divergence corresponds to compression, while positive divergence corresponds to rarefaction. 

Open with DEXTER 
3 The probability density function of the gas density
It is interesting to study the probability distribution of turbulent density fluctuations, because it is a key ingredient for many analytical models of star formation: it is used to explain the stellar initial mass function (Hennebelle & Chabrier 2009,2008; Padoan & Nordlund 2002), the star formation rate (Padoan & Nordlund 2009; Krumholz & McKee 2005; Krumholz et al. 2009), the star formation efficiency (Elmegreen 2008), and the KennicuttSchmidt relation on galactic scales (Elmegreen 2002; Tassis 2007; Kravtsov 2003).
The probability to find a volume with gas density in the range is given by the integral over the volumeweighted probability density function (PDF) of the gas density: . Thus, the PDF p describes a probability density, which has dimensions of probability divided by gas density in the case of . By the same definition, p_{s}(s) denotes the PDF of the logarithmic density .
Figure 4 presents the comparison of the timeaveraged volumeweighted density PDFs p_{s}(s) obtained for solenoidal and compressive forcing. The linear plot of p_{s}(s) (top panel) displays the peak best, whereas the logarithmic representation (bottom panel) reveals the low and highdensity wings of the distributions. Three different fits to analytic expressions (discussed below) are shown as well.
3.1 The density PDF for solenoidal forcing
In numerical experiments of driven supersonic isothermal turbulence
with solenoidal and/or weakly compressive forcing (e.g., Mac Low 1999;
Padoan
et al. 1997; Nordlund & Padoan 1999;
Li
et al. 2003; Beetz et al. 2008; Stone
et al. 1998; Padoan
et al. 2004; VázquezSemadeni 1994; Boldyrev
et al. 2002; Kritsuk et al. 2007),
but also in decaying turbulence (e.g., Ostriker et al. 2001;
Glover
& Mac Low 2007b; Klessen 2000; Ostriker
et al. 1999) it was shown that the density
PDF p_{s}
is close to a lognormal distribution,
where the mean is related to the standard deviation by due to the constraint of mass conservation (e.g., VázquezSemadeni 1994):
Equation (11) simply states that the mean density has to be recovered. This constraint together with the PDF normalisation,
must always be fulfilled for any density PDF whether lognormal or nonGaussian.
Figure 3: Top panel: minimum and maximum logarithmic density as a function of the dynamical time T. Note that compressive forcing yields much stronger compression and rarefaction compared to solenoidal forcing, although the rms Mach number is roughly the same in both cases (see Federrath et al. 2009, Fig. 2). Bottom panel: rms vorticity and rms divergence as a function of the dynamical time. Within the first 2 T, a statistically steady state was reached for both solenoidal (sol) and compressive (comp) forcing. This allows us to average statistical measures (probability density functions, centroid velocity increments, principal component analysis, Fourier spectra and variances) in the range to improve statistical significance of our results and to estimate the amplitude of temporal fluctuations (snapshottosnapshot variations) between different realisations of the turbulence. 

Open with DEXTER 
Figure 4: Volumeweighted density PDFs p(s) of the logarithmic density in linear scaling ( top panel), which displays the peak best, and in logarithmic scaling ( bottom panel) to depict the low and highdensity wings. The PDF obtained from compressive forcing ( ) is significantly wider than the solenoidal one ( ). The peak is shifted to lower values of the logarithmic density s, because of mass conservation, defined in Eq. (11). The density PDF from solenoidal forcing is compatible with a Gaussian distribution. However, there are also nonGaussian features present, which are associated with intermittency effects. These are more prominent in the density PDF obtained from compressive forcing, exhibiting statistically significant deviations from a perfect lognormal (fit using Eq. (10) shown as dashed lines). A skewed lognormal fit (dashdotted lines) given by Eq. (14) provides a better representation, but still does not fit the highdensity tail of the PDF obtained for compressive forcing. Both the PDF data obtained from solenoidal and compressive forcing are best described as lognormal distributions with higherorder corrections defined in Eq. (17), which take into account both the nonGaussian skewness and kurtosis of the distributions. These fits are shown as solid lines (skewkurtlognormal fit). The first four standardised moments defined in Eqs. (13) of the distributions in and s are summarised in Table 1 together with the fit parameters. The grey shaded regions indicate 1 error bars due to temporal fluctuations of the distributions in the regime of fully developed, supersonic turbulence. A total number of 1024^{3} data points contribute to each PDF. 

Open with DEXTER 
From our simulations, we obtain density PDFs in agreement with lognormal distributions for solenoidal forcing. The lognormal fit using Eq. (10) is shown in Fig. 4 as dashed lines. However, the PDF is not perfectly lognormal, i.e., there are weak nonGaussian contributions (see also, Dubinski et al. 1995), especially affecting the wings of the distribution. The strength of these nonGaussian features is quantified by computing higherorder moments (skewness and kurtosis) of the distributions. The first four standardised central moments (see, e.g., Press et al. 1986) of a discrete dataset with N elements are defined as
Note that in our definition of the kurtosis (also called flatness), the Gaussian distribution has . We have computed the first four statistical moments of the volumetric density PDFs shown in Fig. 4. The results are summarised in Table 1. The 1 error given for each statistical moment was obtained by averaging over 81 realisations of the turbulence as described in Sect. 2.2. Both solenoidal and compressive forcing yield density PDFs with deviations from the Gaussian 3rd order (skewness ) and 4th order (kurtosis ) moments.
Table 1: Statistical moments and fit parameters of the PDFs of the volumetric density for solenoidal and compressive forcing shown in Fig. 4.
3.2 The density PDF for compressive forcing
Contrary to the solenoidal case, the PDF obtained for compressive
forcing is not at all well fitted with the perfect lognormal
functional form (dashed line in Fig. 4 for
compressive forcing). Due to the constraints of mass conservation
(Eq. (11))
and normalisation (Eq. (12)),
the peak position and its amplitude cannot be reproduced
simultaneously. The skewness and kurtosis for the compressive forcing
case are also listed in Table 1.
NonGaussian values of skewness and kurtosis, i.e., higherorder
moments require modifications to the analytic expression of the
lognormal PDF given by Eq. (10).
A first step of modification is to allow for a finite
skewness, which is possible with a skewed lognormal distribution (Azzalini 1985)
where , and are fit parameters. Defining , the first four standardised central moments of the distribution are linked to the parameters , and , such that
Skewed lognormal fits are added to Fig. 4 as dashdotted lines and the corresponding fit parameters are given in Table 1. However, for a skewed lognormal distribution, the kurtosis is a function of the skewness, since the skewness and kurtosis in Eqs. (15) both depend on the same parameter only.
Better agreement between an analytic functional form and the
measured PDF can be obtained, if the actual kurtosis of the
data is taken into account as an independent parameter in the
analytical approach. The fundamental derivation of a standard Gaussian
distribution is given by
(16) 
where one parameter is constrained by the normalisation and the two remaining ones are determined by the mean and the dispersion. We can extend this to a modified Gaussianlike distribution by including higherorder moments:
Here, the expansion is stopped at the 4th moment. One parameter is again given by the normalisation, and the remaining four parameters are related to the mean, dispersion, skewness and kurtosis. Fits obtained with this formula are included in Fig. 4 as solid lines. The fit parameters are listed in Table 1. This new functional form is in good agreement with the data from solenoidal and compressive forcing, fitting both the peak and the wings very well. They follow the constraints of mass conservation and normalisation given by Eqs. (11) and (12). We have computed the first four moments of the fitted function and find very good agreement with the first four moments of the actual PDFs.
The fitted parameters a_{3} and a_{4}, which represent the higherorder terms tend to zero compared to the standard Gaussian parameters a_{0}, a_{1} and a_{2} (see Table 1). This means that the higherorder corrections to the standard Gaussian are small. However, we point out that they are absolutely necessary to obtain a good analytic representation of the PDF data, given the fact that Eqs. (11) and (12) must always be fulfilled and that the analytic PDF should return the correct values of the numerically computed moments of the measured distributions.
In the various independent numerical simulations mentioned above, the density PDFs were close to lognormal distributions as in our solenoidal and compressive forcing cases. However, most of these studies also report considerable deviations from Gaussian PDFs, which affected mainly the low and highdensity wings of their distributions. These deviations can be associated with rare events caused by strong intermittent fluctuations during headon collisions of strong shocks and oscillations in very lowdensity rarefaction waves (e.g., Passot & VázquezSemadeni 1998; Kritsuk et al. 2007). The pronounced deviations from the lognormal shape of the density PDF for compressively driven turbulence were also discussed by Schmidt et al. (2009). Even stronger deviations from lognormal PDFs were reported in strongly selfgravitating turbulent systems (e.g., Federrath et al. 2008a; Kainulainen et al. 2009; Klessen 2000).
Intermittency is furthermore inferred from observations, affecting the wings of molecular line profiles (Falgarone & Phillips 1990), and the statistics of centroid velocity increments (HilyBlant et al. 2008). Goodman et al. (2009) measured column density PDFs using dust extinction and emission, as well as molecular lines of gas in the Perseus MC. Using dust extinction maps, Lombardi et al. (2006) obtained the column density PDF for the Pipe nebula. The PDFs found in these studies roughly follow lognormal distributions. However, deviations from perfect lognormal distributions are clearly present in the density PDFs obtained in these studies. They typically exhibit nonGaussian features. For instance, Lombardi et al. (2006) had to apply combinations of multiple Gaussian distributions to obtain good agreement with the measured PDF data.
Figure 5: Volumeweighted correlation PDFs of local Mach number M versus logarithmic density s for solenoidal ( left) and compressive forcing ( right). Adjacent contour levels are spaced by in probability density. Density and Mach number exhibit a very weak, but nonzero correlation in both forcing cases, which provides an explanation for the nonGaussian features in the density PDFs of Fig. 4 (VázquezSemadeni 1994; Passot & VázquezSemadeni 1998; Kritsuk et al. 2007). The two solid lines, intersecting the maxima of both distributions show the mean Mach number as a function of the logarithmic density . The tendency for highdensity gas having lower Mach numbers on average is indicated as power laws in the highdensity parts of the distributions. This suggests that the Mach number for solenoidal and for compressive forcing. 

Open with DEXTER 
3.3 DensityMach number correlation and signatures of intermittency in the density PDFs
As discussed by Passot & VázquezSemadeni (1998), a Gaussian distribution in the logarithm of the density, i.e., a lognormal distribution in is expected for supersonic, isothermal turbulent flows. The fundamental assumption behind this model is that density fluctuations are built up by a hierarchical process. The local density at a given position is determined by a Markov process, i.e., by the product of a large number of independent random fluctuations in time (VázquezSemadeni 1994). If these fluctuations were indeed independent, the quantity would be determined by the sum of this large number of local fluctuations and the distribution in s becomes a Gaussian distribution according to the central limit theorem. Since the Eqs. (2) and (3) are invariant under the transformation for an arbitrary constant s_{0}, the random variable s(t_{n}) should be independent of the local Mach number, and independent of the density at previous times . As pointed out by VázquezSemadeni (1994) and Passot & VázquezSemadeni (1998), this independence breaks down in strong shocks and density extrema, because s_{0} cannot be arbitrarily high due to mass conservation, and an upper boundary s_{+} exists. In consequence, if s_{+} is reached locally, the density cannot increase anymore by a subsequent fluctuation, and the next density is not independent of the previous timestep, causing the fundamental assumption to break down. This also applies to strong rarefaction waves, because creating shocks always produces strongly rarefied regions outside the shock.When the fundamental assumption breaks down, density and velocity statistics are expected to become correlated (VázquezSemadeni 1994; Passot & VázquezSemadeni 1998; Kritsuk et al. 2007). Since in isothermal gas, the sound speed is constant, this translates directly into Mach numberdensity correlations. The average local Mach number M = may therefore exhibit some dependence on the average local density. For instance, it is intuitively clear that headon collisions of strong shocks produce very high density peaks. In the stagnation point of the flow, the local velocity and consequently the local Mach number will almost drop to zero. The time evolution of the maximum and minimum density in Fig. 3 shows these intermittent fluctuations (see also, Porter et al. 1992b; Kritsuk et al. 2007). The intermittent phenomenon corresponds to the situation explained above, for which s_{+} might have been reached, and some dependence of the Mach number on density is expected.
In real molecular clouds, the maximum densities are similarly bounded, and cannot reach infinitely high values, either. This is  unlike the finite resolution constraints in simulations  because the gas becomes optically thick at a certain density ( ), and cannot cool efficiently anymore (e.g., Jappsen et al. 2005; Larson 1969,2005; Penston 1969, and references therein). The gas is not close to isothermal anymore in this regime, and adiabatic compression induced by turbulent motions remain finite in real molecular clouds. Thus, the reason for the breakdown of the densityMach number independence is different in simulations and observations, but it might still be fundamental for the deviations from a lognormal PDF. Moreover, the existence of a characteristic scale may lead to a breakdown of the hierarchical model, and thus to a breakdown of the fundamental assumption. The scale at which supersonic turbulence becomes subsonic is such a scale. This scale is called the sonic scale, and is discussed later in Sect. 8.
We have computed the probability distributions for Mach numberdensity correlations. Figure 5 shows the volumeweighted correlation PDFs of local Mach number M versus density . Although the correlation between density and Mach number is weak as expected for isothermal turbulence (VázquezSemadeni 1994; Passot & VázquezSemadeni 1998), these two quantities are not entirely uncorrelated, which may explain the deviations from perfect lognormal distributions. There is a weak trend for highdensity regions to exhibit lower Mach numbers on average. Powerlaw estimates for densities above the mean logarithmic density indicate Mach numberdensity correlations of the form for solenoidal and for compressive forcing. A similar power law exponent can be obtained from Kritsuk et al. (2007, Fig. 4).
3.4 Numerical resolution dependence of the density PDFs
The highdensity tails of the PDFs in Fig. 4 are not perfectly fit, even when the skewness and kurtosis are taken into account. This is partly due to nonzero 5th, 6th and higherorder moments in the distributions, and partly because our numerical resolution is insufficient to sample the highdensity tail perfectly. Figure 6 shows that even at a numerical resolution of 1024^{3} grid points, the highdensity tails are not converged in both solenoidal and compressive forcing and tend to underestimate high densities. This limitation is shared among all turbulence simulations (see, for instance, the turbulence comparison project by Kitsionas et al. 2009), since the strongest and most intermittent fluctuations building up in the tails will always be truncated due to limited numerical resolution (see also Kowal et al. 2007; Price & Federrath 2010; Hennebelle & Audit 2007). However, the peak and the standard deviation of the PDFs are reproduced quite accurately at a resolution of 256^{3}. Table 2 shows the values of the linear standard deviation and logarithmic standard deviation for numerical resolutions of 256^{3}, 512^{3} and 1024^{3}. There appears to be no strong systematic dependence of the standard deviations on the numerical resolution for resolutions above 256^{3}. The statistical fluctuations are the dominant source of uncertainty in the derived values of the standard deviations. It should be noted however that we have tested only the case of an rms Mach number of about 56 here. There might be a stronger resolution dependence for higher Mach numbers, due to the stronger shocks produced in higher Mach number turbulence, which should be tested in a separate study.
Figure 6: Density PDFs at numerical resolutions of 256^{3}, 512^{3} and 1024^{3} grid cells. The PDFs show very good overall convergence, especially around the peaks. Table 2 shows that the standard deviations are converged with numerical resolution. The highdensity tails, however, are not converged even at a numerical resolution of 1024^{3} grid points, indicating a systematic shift to higher densities with resolution. This limitation is shared among all turbulence simulations (see also, Kitsionas et al. 2009; Price & Federrath 2010; Hennebelle & Audit 2007). The lowdensity wings are subject to strong temporal fluctuations due to intermittent bursts caused by headon collisions of shocks followed by strong rarefaction waves (e.g., Kritsuk et al. 2007). The intermittency causes deviations from a perfect Gaussian distribution and accounts for nonGaussian higherorder moments (skewness and kurtosis) in the distributions. 

Open with DEXTER 
Table 2: Standard deviations of the density PDFs as a function of numerical resolution for solenoidal and compressive forcing shown in Fig. 6.
3.5 The column density PDFs and comparison with observations
The strong difference between the statistics of the solenoidal and compressive forcing cases seen in the PDFs of the volumetric density shown in Fig. 4 is reflected by the corresponding column density PDFs. The timeaveraged and projectionaveraged column density PDFs are shown in Fig. 7. Analogous to Table 1 for the volumetric density PDFs, we summarise the statistical quantities and fit parameters for the column density PDFs in Table 3. The main results and conclusions obtained for the volumetric density distributions also hold for the column density distributions. Compressive forcing yields a column density standard deviation roughly three times larger than solenoidal forcing. The relative difference between solenoidal and compressive forcing is thus roughly the same for the volumetric and the column density distributions. However, the absolute values are lower for the column density distributions compared to the volumetric density distributions. The reason for this is that by computing projections of the volumetric density fields, density fluctuations are effectively averaged out by integration along the lineofsight, and as a consequence, the column density dispersions become smaller compared to the corresponding volumetric density dispersions.
Table 3: Same as Table 1, but for the PDFs of the column density shown in Fig. 7.
The small inset in the upper right corner of Fig. 7 additionally shows the column density PDFs computed along the zaxis at one single time t=2 T corresponding to the map shown in Fig. 2. This figure shows the effect of studying one realisation only, without time and/or projectionaveraging. This is interesting to consider, because observations can only measure column density distributions at one single time. Improving the statistical significance would only be possible by studying multiple fields and averaging in space rather than in time invoking the ergodic theorem as suggested by Goodman et al. (2009). However, even by studying one turbulent realisation only, the difference between solenoidal and compressive forcing is recovered from the dispersions of the distributions. We therefore expect that using observations of column density PDFs, one can distinguish purely solenoidal from purely compressive forcing by measuring the dispersion of the column density PDF.
Figure 7: Same as Fig. 4, but the time and projectionaveraged logarithmic column density PDFs of are shown. and denote the column density and the mean column density integrated along the i=x, y, z principal axes respectively. As for the volumetric PDFs of Fig. 4, the standard deviation of the column density PDF obtained from compressive forcing is roughly three times larger than from solenoidal forcing (see Table 3). The inset in the upper right corner shows the PDFs of column density computed in zprojection at a fixed time t=2 T, corresponding to the snapshots shown in Fig. 2. The density dispersions computed for these instantaneous PDFs are and for solenoidal forcing, and and for compressive forcing. Although these distributions are quite noisy, the influence of the forcing is still clearly discernible. Thus, by studying instantaneous column density PDFs, which are accessible to observations, one should be able to distinguish solenoidal from compressive forcing. 

Open with DEXTER 
Figure 8: Diamonds: the proportionality parameter b in the density dispersionMach number relation, Eq. (18), computed as for eleven 3D models at numerical resolution of 256^{3} grid cells ( top panel) and eleven 2D models at numerical resolution of 1024^{2} grid cells ( bottom panel), ranging from purely compressive forcing () to purely solenoidal forcing (). The parameter b decreases smoothly from for compressive forcing to in 3D and in 2D for solenoidal forcing. Stars: ratio of longitudinal to total power in the velocity power spectrum (see Sect. 6.1). This quantity provides a measure for the relative amount of compression induced by the turbulent velocity field, and appears to be correlated with the standard deviation of the density PDF. Squares: same as stars, but multiplied by the geometrical factor with D=3 for the threedimensional case and D=2 for the twodimensional case. The quantity provides a good numerical estimate of the PDF parameter b. The dashed lines show model fits using Eq. (23) for D=3 ( top panel) and D=2 ( bottom panel). 

Open with DEXTER 
Goodman et al. (2009) provided measurements of the column density PDFs in the Perseus MC obtained with three different methods: dust extinction, dust emission, and ^{13}CO gas emission. Although systematic differences were found between the three methods, they conclude that in general, the measured column density PDFs are close to, but not perfect lognormal distributions, which is consistent with our results. They furthermore provided the column density PDFs and the column density dispersions for six subregions in the Perseus MC. The difference between the dispersions measured for these subregions is not as large as the difference between purely solenoidal and purely compressive forcing. The largest difference in the column density dispersions among the six subregions found by Goodman et al. (2009) is only about 50% relative to the average column density dispersion measured in the Perseus MC. This indicates that both purely solenoidal and purely compressive forcing are very unlikely to occur in nature. On the other hand, a varying mixture of solenoidal and compressive modes close to the natural mixture of 2:1 can easily explain the 50% difference in density dispersion measured among the different regions. In particular, the Shell region (Ridge et al. 2006), which surrounds the B star HD 278942 exhibits the largest density dispersion among all the subregions studied by Goodman et al. (2009), although its velocity dispersion is rather small compared to the others. This indicates that turbulent motions may be driven compressively rather than solenoidally within the Shell region. Goodman et al. (2009) indeed mentioned that the gas in the Shell is dominated by an ``obvious driver'', skewing the column density distribution towards lower values compared to the other regions. Due to the constraints of mass conservation (Eq. (11)) and normalisation (Eq. (12)), both the peak position and the peak value of the PDF skew to lower values, if the density dispersion increases (see Fig. 7). Taken together, this suggests that the Shell in the Perseus MC represents an example of strongly compressive turbulence forcing rather than purely solenoidal forcing.
3.6 The forcing dependence of the density dispersionMach number relation
In Federrath et al. (2008b), we investigated the density dispersionMach number relation (Padoan et al. 1997; Passot & VázquezSemadeni 1998)^{},This relation was also investigated in Kowal et al. (2007, Fig. 11), indicating that the standard deviation of turbulent density fluctuations, is directly proportional to the sonic Mach number in the supersonic regime. It must be noted, however, that it was only directly tested for rather low rms Mach numbers, (Kowal et al. 2007) and (Passot & VázquezSemadeni 1998), compared to typical Mach numbers in molecular clouds. If additionally a lognormal PDF, Eq. (10) is assumed, then Eq. (18) can be expressed as
with the same parameter b (Federrath et al. 2008b; Padoan et al. 1997).
We begin our discussion of the forcing dependence of the density dispersionMach number relation with a problem raised by Mac Low et al. (2005) and Glover & Mac Low (2007b). Mac Low et al. (2005) and Glover & Mac Low (2007b) claimed that the density dispersionMach number relation found by Passot & VázquezSemadeni (1998), (which is a Taylor expansion of Eq. (19) for small rms Mach numbers), with did not at all fit their results for pressure and density PDFs, while Eq. (19) with (Padoan et al. 1997) provided a much better representation of their data. The main difference in the density dispersionMach number relations by Padoan et al. (1997) and Passot & VázquezSemadeni (1998) is the proportionality constant b. It is and in Padoan et al. (1997) and Passot & VázquezSemadeni (1998), respectively. Our forcing analysis provides the solution to this apparent difference, which lies at the heart of the disagreement of the PDF data analysed in Mac Low et al. (2005) and Glover & Mac Low (2007b) with the model by Passot & VázquezSemadeni (1998). Passot & VázquezSemadeni (1998) used 1D models. In 1D, only compressive forcing is possible, because no transverse waves can exist. In contrast, Mac Low et al. (2005) and Glover & Mac Low (2007b) used a mixture of solenoidal and compressive forcing in 3D. In this section, we show that the parameter b in both Eqs. (18) and (19) is a function of the forcing parameter . Indeed, using the relation analysed in Passot & VázquezSemadeni (1998), but with a lower proportionality constant (b=0.5 in contrast to b=1) gives a very good representation of the PDF data in Mac Low et al. (2005, Fig. 8). Thus, an investigation of the parameters that control b seems necessary and important.
Moreover, relations (18)
and (19)
are key ingredients for the analytical models of the stellar initial
mass function by Padoan
& Nordlund (2002) and Hennebelle & Chabrier
(2008), as well as for the star formation rate model
by Krumholz & McKee
(2005) and Krumholz
et al. (2009) and for the star formation efficiency
model by Elmegreen (2008).
In all these models, b is assumed to
be 0.5, which is an empirical result from
magnetohydrodynamical simulations by Padoan et al. (1997).
On the other hand, Passot
& VázquezSemadeni (1998) found
from 1D hydrodynamical simulations. Federrath
et al. (2008b) resolved this disagreement between Padoan et al. (1997)
and Passot &
VázquezSemadeni (1998) by showing that b
is a function of the ratio
of compressive to solenoidal
modes of the turbulence forcing. However, Federrath
et al. (2008b) only tested the two limiting cases of
purely solenoidal forcing ()
and purely compressive forcing (). They approximated the
regime of mixtures with a heuristic model, which had a linear
dependence of b on :
Here, we refine this model based on eleven additional simulations with separated by for rms Mach numbers of 5 in 2D and 3D. These simulations allow us to test eleven different mixtures of forcing controlled by the parameter (see Eq. (9)). The models were run at a numerical resolution of 256^{3} and 1024^{2} grid points in 3D and 2D, respectively. We use a lower resolution in 3D, because using our standard resolution of 1024^{3} would be too computationally intensive. However, as shown in Sect. 3.4, the standard deviation of the density is fairly well reproduced at 256^{3}, as is the rms Mach number (Federrath et al. 2009, Fig. 2), which allows a reasonably accurate determination of b. The results are plotted in Fig. 8 as diamonds for 3D (top panel) and 2D (bottom panel). We used Eq. (18) to measure b, because unlike Eq. (19), this version of the standard deviationMach number relation does not rest on the assumption of a lognormal PDF. In fact, if Eq. (19) was used to derive b for models with , b would be overestimated significantly (by up to an order of magnitude for ), because the deviations from the perfect lognormal distribution are stronger for (cf. Sect. 3.2; see also Schmidt et al. (2009)).
Figure 8 shows that the dependence of b on is nonlinear. For 3D turbulence the parameter b increases smoothly from for to for , and for 2D turbulence from for to for . However, there is an apparent break at , which represents the natural forcing mixture used in many previous studies. For the bparameter remains close to the value obtained for purely solenoidal forcing, i.e. in 3D and in 2D. The flat part of the data in Fig. 8 for explains why in previous studies with a natural forcing mixture (e.g., Klessen et al. 2000; Glover et al. 2010; Li et al. 2003; Mac Low et al. 1998; Kritsuk et al. 2007), the turbulence statistics were close to the purely solenoidal forcing case (e.g., Kowal et al. 2007; Stone et al. 1998; Padoan et al. 1997; Padoan & Nordlund 2002; Boldyrev et al. 2002; Burkhart et al. 2009; Lemaster & Stone 2008). In contrast, b increases much more strongly for , until it reaches for purely compressive forcing (e.g., Federrath et al. 2008b; Schmidt et al. 2009; Passot & VázquezSemadeni 1998).
Equation (20)
thus needs to be refined to account for the nonlinear dependence
of b on the forcing. Moreover,
Eq. (20)
was based on the analytic expression of the forcing parameter
(cf. Sect. 2.1).
However, the numerical estimate of b
depends on how well the code can actually induce compression through
the buildup of divergence in the velocity field. Thus, different codes
can produce slightly different values of b
for the same forcing parameter .
This is because of the varying efficiency of codes to convert the
energy provided by a given forcing into actual velocity fluctuations
(e.g. Kitsionas
et al. 2009; Price & Federrath 2010).
To construct a refined model for b
that does not directly rest on the analytic forcing parameter
and that accounts for the nonlinear dependence on the forcing, we
recall that b is a normalised measure of
compression. Compression is caused by converging flows and shocks,
which have a finite magnitude of velocity divergence.
A normalised measure of compression is thus also provided by
dividing the power in longitudinal modes of the velocity field by the
total power of all modes in the velocity field,
(21) 
We therefore expect a dependence of b on .
Figure 8
shows
as a function of
(plotted as stars) for 3D and 2D turbulence.
It is indeed correlated with b,
however,
is less than b by a factor of
roughly
in 3D and
in 2D. The squares in Fig. 8 show
in 3D
and
in 2D, which seems to provide a good estimate of b.
The factor
is a geometrical factor for 3D turbulence
(the diagonal in a cube of size unity).
It is
in 2D turbulence (the diagonal in a square of size
unity), and
in 1D. The latter in particular is trivial, because
in 1D only longitudinal modes can exist, and thus
for any value of
(cf. Fig. 1). The
larger geometrical factors in 2D and 3D account for
the fact that the longitudinal velocity fluctuations, which induce
compression occupy only one of the available spatial directions
(two in 2D and three in 3D) on average. For
the general case of supersonic turbulence in D=1,2
and 3 dimensions, these ideas lead to
which is solely based on the ratio of the power in longitudinal modes in the velocity field to the total power of all modes in the velocity field, .
In addition to the refined model based on the compressive
ratio
in
Eq. (22),
we provide a fit function for b based on
the forcing parameter .
The dashed lines in Fig. 8 show
The forcing ratio is given by Eq. (9). The first summand in Eq. (23) is the expected ratio of longitudinal modes (compression) in a supersonic turbulent medium for a purely solenoidal forcing, i.e. a forcing that does not directly induce compression. The second summand is the contribution to the compression directly induced by the forcing. The model Eq. (23) is similar to Eq. (20), but with a nonlinear dependence of b on the forcing parameter .
We suggest that the dependence of b on the forcing solves a puzzle reported by Pineda et al. (2008). They provided measurements of velocity dispersions and ^{12}CO excitation temperatures for the six subregions in the Perseus MC. The molecular excitation temperatures serve as a guide for the actual gas temperature, from which the sound speed can be estimated. From these values, the local rms Mach numbers are computed as the ratio of the local velocity dispersion to the local sound speed. Goodman et al. (2009) and Pineda et al. (2008) pointed out that there is clearly no correlation of the form suggested by Eq. (19) for a fixed parameter b across the investigated subregions in the Perseus MC. For instance, the Shell region exhibits an intermediate to small velocity dispersion derived from ^{12}CO and ^{13}CO observations, while its density dispersion is the largest in the Perseus MC. This provides additional support to our suggestion that the Shell in Perseus is dominated by compressive turbulence forcing for which b takes a higher value compared to solenoidal forcing. The apparent lack of density dispersionMach number correlation reported by Pineda et al. (2008) and Goodman et al. (2009) for a fixed parameter b can thus be explained, because b is in fact not fixed across different subregions in the Perseus MC.
We plan to measure b in different regions of the ISM in future studies. However, the main problem in a quantitative analysis of Eq. (18) with observational data is that the column density dispersion is typically smaller than the 3D density dispersion (compare Tables 1 and 3). The relation between the column density PDF and the volumetric density PDF is nontrivial and depends on whether the column density tracer is optically thin or optically thick and on the scale of the turbulence driving. However, Brunt et al. (2010) developed a promising technique to estimate the 3D density variance from 2D observations with an accuracy of about 10%.
4 Intermittency
Intermittency manifests itself in
 i)
 nonGaussian (often exponential) wings of PDFs of quantities involving density and/or velocity, its derivatives (e.g., vorticity) and combinations of density and velocity (e.g., and as discussed in Appendix A);
 ii)
 anomalous scaling of the higherorder structure functions of the velocity field (e.g., Anselmet et al. 1984) and centroid velocity increments (HilyBlant et al. 2008; Lis et al. 1996); and
 iii)
 coherent structures of intense vorticity ( ) (see Moisy & Jiménez 2004; Vincent & Meneguzzi 1991, for results of incompressible turbulence), and of strong shocks and rarefaction waves ( ).
Figure 9: PDFs of centroid velocity increments, computed using Eqs. (24) and (25) are shown as a function of lag in units of grid cells for solenoidal forcing ( left) and compressive forcing ( right). The PDFs are very close to Gaussian distributions for long lags, whereas for short lags, they develop exponential tails, which is a manifestation of intermittency (e.g., HilyBlant et al. 2008, and references therein). 

Open with DEXTER 
4.1 The probability distribution of centroid velocity increments
Since there is evidence of filamentary coherent structures in the
vorticity (intermittency item iii) of our models, and because
there is additional evidence of nonGaussian tails in the density PDFs
(intermittency item i) discussed in Sect. 3, we now proceed
to examine the PDFs and the scaling of centroid velocity increments
(intermittency item ii) to assess the strength of the
intermittency. We compare centroid velocity increments (CVIs)
for solenoidal and compressive forcing and discuss the interpretation
of observations based on that comparison. Following the analysis by Lis et al. (1996),
who discuss CVIs computed for the turbulence simulation by Porter et al.
(1994), and following the CVI analysis of the
Polaris Flare and of the Taurus MC by HilyBlant
et al. (2008), the centroid velocity increment is
defined as
where the angle average is computed over all possible directions of the vector in the plane perpendicular to the lineofsight. Thus, only depends on the norm of the lag vector , which separates two points and in the plane of the sky (x,y). The normalised centroid velocity, in Eq. (24) is defined as
The variable denotes the lineofsight velocity in zdirection. We have however computed separately along each of the three principal linesofsight x, y and z of our Cartesian domain in order to examine the effects of varying the projection. Also note that we have computed normalised centroid velocities (Lazarian & Esquivel 2003), since we want to compare to HilyBlant et al. (2008). Another point to mention here is that the centroid velocities, are typically computed using an intensity weighting instead of a density weighting. This is because the gas density cannot be measured directly, whereas the emission intensity is accessible to observations. By using density weighting we implicitly assume optically thin emission. For optically thick emission, uniform weighting would be more appropriate (Lis et al. 1996).
Figure 9 shows the PDFs of computed for varying lag in units of the numerical cell size . They should be compared to HilyBlant et al. (2008, Fig. 46). The PDFs are mainly Gaussian for large lags, whereas for smaller separations, they develop exponential tails, indicating intermittent behaviour. This result is consistent with the numerical simulation analysed by Lis et al. (1996), and with observations of the Oph Cloud, the Orion B and the Polaris Flare by Lis et al. (1998), Miesch et al. (1999) and HilyBlant et al. (2008), respectively.
Following the analysis by HilyBlant et al. (2008), we computed the kurtosis of the PDFs of CVIs using the definition in Eqs. (13). Note that corresponds to a Gaussian distribution, and corresponds to an exponential function. The kurtosis of the CVI PDFs is shown in Fig. 10 as a function of spatial lag , and can be directly compared to HilyBlant et al. (2008, Fig. 7). Both forcing types exhibit nearly Gaussian values of the kurtosis at lags . On the other hand, for , both forcing types produce nonGaussian PDFs. Solenoidal forcing approaches the exponential value for . Compressive forcing yields exponential values already for lags , while solenoidal forcing has on these scales. This indicates stronger intermittency in the case of compressive forcing. For , compressive forcing yields even superexponential values of . For both solenoidal and compressive forcings, we show later in Sect. 6 that is in the dissipation range for numerical turbulence. Compressive velocity modes dominate in this regime (see Fig. 14), which may result artificially in extreme intermittency. For , compressive forcing gives 1.0, which is roughly 35% larger than the Polaris Flare observations at their resolution limit. The solenoidal case on the other hand gives 0.5, which is in very good agreement with the IRAM and KOSMA data discussed by HilyBlant et al. (2008, Fig. 7). Depending on the actual lag used for the comparison, both solenoidal and compressive forcing seem to be consistent with the observations. However, it should be noted that the lags cannot be easily compared for the real clouds and the simulations, because simulated and observed fields have different spatial resolution. Moreover, the simulated fields have periodic boundaries, while the true fields don't. Nevertheless, the similarity of the observed and the numerically simulated CVIs indicates that turbulence intermittency plays an important role in both our simulations and in real molecular clouds.
Figure 10: Kurtosis of the PDFs of centroid velocity increments shown in Fig. 9 as a function of the lag in units of grid cells for solenoidal and compressive forcing. Note that a kurtosis value of 3 (horizontal dotdashed line) corresponds to the value for a Gaussian distribution. NonGaussian values of the kurtosis are obtained for . The error bars contain both snapshottosnapshot variations as well as the variations between centroid velocity increments computed by integration along the x, y and z axes. This figure can be compared to observations of the Polaris Flare and Taurus MC (see Fig. 7 of HilyBlant et al. 2008). 

Open with DEXTER 
Figure 11: Scaling of the structure functions of centroid velocity increments defined in Eq. (26) for solenoidal forcing ( left) and compressive forcing ( right) up to the 6th order. Scaling exponents obtained using powerlaw fits following Eq. (27) within the inertial range are indicated in the figures and summarised in Table 4. 

Open with DEXTER 
Table 4: Scaling of the structure functions of centroid velocity increments.
The Polaris Flare has a very low star formation rate and is therefore appropriate for studying the statistics of interstellar supersonic turbulence without contamination by internal energy sources. In contrast, the Taurus MC is actively forming stars. Against our expectations, the Taurus MC data display very weak intermittent behaviour and the kurtosis remains at the Gaussian values in HilyBlant et al. (2008, Fig. 7). However, the Taurus field studied by HilyBlant et al. (2008) is located far from starforming regions in a translucent part of the Taurus MC (Falgarone 2009, private communication). This may explain why the Taurus field displays only very weak intermittency. It would be interesting to repeat the analysis of centroid velocity increments for regions of confirmed star formation, including regions with winds, outflows and ionisation feedback from young stellar objects to see whether these regions indeed display stronger intermittency.
4.2 The structure function scaling of centroid velocity increments
In this section, we discuss the scaling of the pth order
structure function of CVIs, defined as
We have averaged over a large enough sample of independent increments that increasing the sample size produced no change in the value of for , which is demonstrated in Appendix B. Figure 11 shows the CVI structure functions for solenoidal and compressive forcing. The CVI structure functions were fit to power laws of the form
within the inertial range^{}, defined equivalently to the study in Federrath et al. (2009). The value for each powerlaw exponent is indicated in Fig. 11 and summarised in Table 4.
For a direct comparison of CVI structure functions with the study by HilyBlant et al. (2008, Fig. 8), we apply the extended selfsimilarity (ESS) hypothesis (Benzi et al. 1993), which states that the inertial range scaling may be extended beyond the inertial range, such that powerlaw fits can be applied over a larger dynamic range. The ESS hypothesis is used by plotting the pth order against the 3rd order (Benzi et al. 1993). These plots are shown in Fig. 12. Indeed, the scaling range is drastically increased using ESS. All ESS data points are consistent with a single power law for each CVI structure function order . We summarise the scaling exponents with and without using the ESS hypothesis in Table 4.
Figure 12: Same as Fig. 11, but using the extended selfsimilarity hypothesis (Benzi et al. 1993), allowing for a direct comparison of the scaling exponents of centroid velocity increments with the study by HilyBlant et al. (2008) for the Polaris Flare and Taurus MC (see Table 4). 

Open with DEXTER 
Table 4
furthermore provides the ESS scaling exponents obtained for the
Polaris Flare (HilyBlant
et al. 2008, Table 3),
as well as the scaling exponents obtained from intermittency
models of the structure function scaling exponents
by She & Leveque (1994) ( ) and Boldyrev (2002) ( ). In these models, the fractal codimension is related to the fractal dimension of the most intermittent structures by . The She & Leveque (1994) model assumes 1D vortex filaments as the most intermittent structures ( ), whereas the Boldyrev (2002) model assumes sheetlike structures with .
For solenoidal forcing, the scaling of the CVI structure
functions using ESS is very similar to the She
& Leveque (1994) model. This model is appropriate for
incompressible turbulence, for which the most
intermittent structures are expected to be filaments (She & Leveque 1994,
).
Interestingly, their model seems to be consistent with the measurements
in the Polaris Flare by HilyBlant
et al. (2008) and with our solenoidal forcing case.
In contrast, the scaling exponents derived for
compressive forcing are better consistent with the intermittency model
by Boldyrev (2002,
Nevertheless, a direct comparison of CVI structure function scaling obtained in numerical experiments and observations can provide useful information to distinguish between different parameters of the turbulence, as for instance different turbulence forcings.
Figure 13: Principal component analysis (PCA) for solenoidal ( left) and compressive forcing ( right). The PCA slopes obtained for solenoidal and compressive forcings are summarised and compared with observations by Heyer et al. (2006) in Table 5. The error bars contain the contribution from temporal variations and from three different projections along the x, y and zaxes. The data were resampled from 1024^{3} to 256^{3} grid points prior to PCA. The resampling speeds up the PCA and has virtually no effect on the inertial range scaling (see e.g., Federrath et al. 2009; Padoan et al. 2006). 

Open with DEXTER 
Table 5: Comparison of measured PCA scaling slopes.
5 Principal component analysis
Principal component analysis (PCA) is a multivariate tool (Murtagh & Heck 1987)
introduced by Heyer
& Schloerb (1997) for measuring the scaling of
interstellar turbulence. It has been used for studying the
structure and scaling in several molecular cloud regions, simulations
and synthetic images (Heyer et al. 2006;
Brunt
& Heyer 2002b; Heyer & Brunt 2004;
Brunt
et al. 2003; Brunt & Heyer 2002a).
PCA can be used to characterise structure on different scales. For best
comparison with observations, we choose to work in
positionpositionvelocity (PPV) space. Since our simulation data are
typically stored in positionpositionposition (PPP) space, we
transformed our PPP cubes into PPV space prior
to PCA. As for the CVIs discussed in the previous
section, we use the approximation of optically thin radiative transfer
to derive radiation intensity. This means that we essentially assume
that the emission is proportional to the gas density.
The PPV data therefore represent a simulated measured
intensity
at spatial position
and spectral position v_{z,j}.
The indices i and j
thus represent the spatial and spectral coordinates respectively.
A detailed description of the PCA technique is given
by Heyer & Schloerb
(1997) and Brunt
& Heyer (2002a). The most important steps necessary
to derive the characteristic length scales and corresponding velocity
scales using PCA are described below. First, the covariance matrix
(29) 
is constructed by summation over all spatial points N_{x} and N_{y}. Solving the eigenvalue equation
(30) 
yields the lth eigenvalue and the lth eigenvector of the covariance matrix. The subsequent projection
(31) 
onto the eigenvectors yields the lth eigenimage I_{i}^{(l)}. Autocorrelation functions (ACFs) are then computed for each of the eigenimages and eigenvectors. The spatial scale on which the twodimensional ACF of the lth eigenimage falls off by 1/e defines the lth characteristic spatial scale. Following the same procedure, the corresponding characteristic velocity scale is determined from the ACF of the lth eigenvector, which contains the spectral information.
Figure 13 shows our time and projectionaveraged set of spatial and velocity scales obtained with PCA. We have fitted power laws to the PCA data, which yielded PCA scaling exponents for solenoidal and compressive forcing respectively. For solenoidal forcing we find 0.05 and for compressive forcing we find 0.09 (see Table 5). The different PCA slopes derived for solenoidal and compressive forcing suggest that by using PCA, differences in the mixture of transverse and longitudinal modes of the velocity field can be detected. However, the difference between solenoidal and compressive forcing is only at the 1 level.
Figure 14: Top panels: total, transverse (rotational) and longitudinal (compressible) velocity Fourier spectra E(k) defined in Eq. (32) and compensated by k^{2} for solenoidal ( left) and compressive forcing ( right). Error bars indicate temporal variations, which account for an uncertainty of roughly of all scaling slopes reported for the inertial range . The inferred inertial range scaling exponents for both solenoidal and compressive forcing are consistent with independent numerical simulations and with observations of the sizelinewidth relation (see text). Note that the transverse part, falls off more steeply than the longitudinal part, for both forcing types in the inertial range. Bottom panels: ratio of the energy in longitudinal velocity modes to the total energy in velocity modes . For solenoidal forcing, we obtain 1/3 in the inertial range (horizontal dashdotted line), because compression can only occur in one of the three spatial dimensions on average (Federrath et al. 2008b; Elmegreen & Scalo 2004). For compressive forcing, this ratio is roughly 1/2, which corresponds to an equipartition of longitudinal and transverse velocity modes. Note however that compressive forcing can compress the gas in all three spatial dimensions directly, whereas solenoidal forcing can only induce compression indirectly through the velocity field (Federrath et al. 2008b). The excess of longitudinal modes at high wavenumbers stems from numerical dissipation, which is more effectively dissipating transverse than longitudinal modes on small scales due to the discretisation onto a grid. This suggests that roughly 30 grid cells are needed to accurately resolve a vortex, while a shock is typically resolved with roughly 3 grid cells using the piecewise parabolic method (Colella & Woodward 1984). However, for a numerical resolution of 1024^{3} grid cells, we find that wavenumbers are almost unaffected by the discretisation and by the parameters of the numerical scheme (see Appendix C). 

Open with DEXTER 
Heyer et al. (2006) applied PCA to the Rosette MC and to G2162.5 (Maddalenas's Cloud). These two clouds are quite different in dynamical and evolutionary state, although they exhibit roughly the same turbulence Mach number. Heyer et al. (2006) measured the Mach number 45 on a scale of for both clouds. The Rosette MC exhibits confirmed massive star formation, whereas G2162.5 has a low star formation rate, similar to the Polaris Flare discussed in the previous section. Heyer et al. (2006) measured PCA slopes for both clouds and additionally provided the PCA slopes in two distinct subregions of the Rosette MC. The first subregion is inside the HII region (Zone I) surrounding the massive star cluster NGC 2244^{}, while the other subregion is outside of this HII region (Zone II). The measured PCA slopes obtained from ^{12}CO and ^{13}CO observations are summarised in Table 5 together with our estimates for solenoidal and compressive forcing. The PCA scaling exponent for solenoidal forcing is very close to the PCA scaling exponents derived from the ^{12}CO observations in the G2162.5 ( 0.04) and in Zone II of the Rosette MC ( 0.06). In contrast, the PCA slope derived from ^{12}CO observations in Zone I of the Rosette MC ( 0.06) is better consistent with our compressive forcing case. This indicates that Zone I contains more kinetic energy in compressive modes than Zone II and G2162.5. The corresponding ^{13}CO observations reported in Heyer et al. (2006) yield slightly larger differences between the PCA scaling exponents derived for Zone I on the one hand, and Zone II and G2162.5 on the other hand (see also Table 5). This supports the idea that Zone I in the Rosette MC, and Zone II as well as G2162.5 contain quite different amounts of compressive modes in the velocity field, which may be the result of different turbulence forcing mechanisms, similar to the differences obtained in purely solenoidal and compressive forcings.
6 Fourier spectra
6.1 Velocity Fourier spectra
Fourier spectra of the velocity field E(k) are typically used to distinguish between Kolmogorov (1941) turbulence, and Burgers (1948) turbulence, . For highly compressible, isothermal, supersonic, turbulent flow, it has been shown that the inertial range scaling is close to Burgers turbulence. For instance, Kritsuk et al. (2007) found and Schmidt et al. (2009) obtained from highresolution numerical simulations.
The Fourier spectrum of a quantity provides a measure of the
scale dependence of this quantity. Velocity Fourier spectra are thus
defined as
where denotes the Fourier transform of the velocity field (e.g., Frisch 1995). The total Fourier spectrum can be separated into transverse ( ) and longitudinal ( ) parts by applying a Helmholtz decomposition. Note that integrating the transverse energy spectrum yields the kinetic energy in transverse (rotational) modes, while integration of the longitudinal energy spectrum yields the kinetic energy in longitudinal (compressible) modes. Furthermore, by integrating the velocity spectrum from k_{1} to k_{2}, one obtains the kinetic energy content on length scales corresponding to the wavenumber interval [k_{1},k_{2}]. Since the mean velocity is zero in our simulations, integration of the total velocity Fourier spectrum E(k) over all wavenumbers yields the total variance of velocity fluctuations :
(33) 
The upper bound of the integral is the cutoff wavenumber for a cubic dataset with N^{3} data points. Thus, for our standard resolution of 1024^{3} grid cells.
In Fig. 14 we show the total velocity Fourier spectra E(k) as defined in Eq. (32) together with its decomposition into transverse and longitudinal parts for solenoidal and compressive forcing respectively. The prominent signature of the different forcings on the main driving scale, k=2 is clearly noticeable: solenoidal forcing excites mostly transverse modes, whereas compressive forcing excites mostly longitudinal modes in the velocity field at k=2. However, the forcing has direct influence only for 1<k<3 (see Sect. 2.1). Further down the cascade, the turbulent flow develops its own statistics as a result of nonlinear interactions in the inertial range . We emphasise that this scaling range was chosen very carefully, since turbulence simulations will only provide a small inertial range even at resolutions of 1024^{3} grid cells (see, e.g., Klein et al. 2007; Lemaster & Stone 2009). This is mainly caused by the bottleneck phenomenon (e.g., Schmidt et al. 2006; Haugen & Brandenburg 2004; Porter et al. 1994; Dobler et al. 2003; Kritsuk et al. 2007), which may slightly affect the Fourier spectra in the dissipation range. However, the bottleneck phenomenon had no significant impact on the turbulence statistics in our numerical study for wavenumbers . This is demonstrated in Appendix C, where we present the resolution dependence of the Fourier spectra and the dependence on parameters of the PPM numerical scheme. We conclude that the statistical quantities derived for wavenumbers are not significantly affected by the numerical scheme or limited resolution applied in the present study.
We apply powerlaw fits to the inertial range data with the resulting slopes indicated in Fig. 14 (top panels). Both solenoidal and compressive forcing yield slopes consistent with sizelinewidth relations inferred from observations (e.g., Myers 1983; Solomon et al. 1987; Perault et al. 1986; Ossenkopf & Mac Low 2002; Heyer et al. 2009; Falgarone et al. 1992; Heyer & Brunt 2004; Padoan et al. 2006; Ossenkopf et al. 2008b; Miesch & Bally 1994; Larson 1981; Padoan et al. 2003), and with the results of independent numerical simulations (e.g., Schmidt et al. 2009; Padoan et al. 2004; Klessen et al. 2000; Boldyrev et al. 2002; Kritsuk et al. 2007). Note that sizelinewidth relations of the form with scaling exponents correspond to Fourier spectra with scaling exponents in the range , because . However, it must be emphasised that the relation between scaling exponents obtained from observational maps of centroid velocities (as discussed in Sect. 4.2) and 3D velocity fields from simulations is nontrivial, because of projectionsmoothing and intensityweighting. Projectionsmoothing increases the scaling exponents of the 2D projection of a 3D field such that (e.g., Brunt & Mac Low 2004; Stutzki et al. 1998). However, Brunt & Mac Low (2004) showed that the effect of projectionsmoothing is compensated statistically (but not identically) by intensityweighting of observed centroid velocity maps. Thus, our measurements of velocity scaling seem consistent with observations.
It is important to note that the transverse parts fall off more steeply than the longitudinal parts for both forcing types within the inertial range. For solenoidal forcing, we find and , and for compressive forcing, and . This result indicates that longitudinal modes can survive down to small scales, such that compression may not be neglected anywhere in the turbulent cascade. Lemaster & Stone (2009, Figs. 9, 10) obtain and for their hydrodynamical model with solenoidal forcing at a resolution of 1024^{3} grid points in the Athena code. This is consistent with our findings for the scale dependence of the transverse and longitudinal parts and shows that the kinetic energy in longitudinal modes must not be neglected within the inertial range.
Figure 15: Left panel: fourier spectra of the velocity, E(k) defined in Eq. (32) (crosses and diamonds) and Fourier spectra of the logarithmic density fluctuations, S(k) defined in Eq. (35) (triangles and squares) for solenoidal and compressive forcing, respectively. Both E(k) and S(k) are compensated by k^{2} allowing for a better determination of the inertial range scaling. The density fluctuation power spectra differ significantly in the inertial range with for solenoidal and for compressive forcing. The scale on which the density fluctuation spectra from solenoidal and compressive forcing cross each other and where the slope obtained in compressive forcing breaks and approaches the shallower slope of the solenoidal forcing case roughly coincides with the sonic wavenumber (vertical dashed lines) defined in Eq. (38). Right panel: same as left panel, but instead of using Fourier spectra to determine the inertial range scaling, we use the variance method to derive the scaling slopes in physical space. Note that the scaling slopes obtained with the variance technique are related to the slopes of the Fourier spectra by (Stutzki et al. 1998). Error bars denote 1 temporal fluctuations. 

Open with DEXTER 
In order to quantify the relative importance of compression over
rotation in the turbulent motions, we present plots of
the ratio
in the bottom panels of Fig. 14. Solenoidal forcing yields 1/3 in the inertial range. We emphasise that the ratio 1/3 was expected from the fact that compression can only occur in one of the three available spatial dimensions on average in the case of supersonic flow driven by a purely solenoidal force field (Federrath et al. 2008b; Elmegreen & Scalo 2004). This is the fundamental idea on which the heuristic model of the density dispersionMach number relation given by Eq. (20) was based. For compressive forcing, we find 1/2 in the inertial range as a result of the direct compression induced by compressive forcing. Thus, solenoidal and compressive forcing produce quite similar amounts of compressive modes in the velocity field ( 1/3 versus 1/2). This means that even fully compressive forcing can behave very similar to solenoidal forcing in the inertial range, as far as pure velocity statistics are concerned. However, we show in the next section that the density statistics are very different in the inertial range. The same is true for combined densityvelocity statistics (see Appendix A).
We also note here that the rise of at for both forcing types is a numerical effect, which comes from the discretisation of the velocity field onto a grid with finite resolution. This shows that energy in rotational modes cannot be accounted for accurately if vortices are smaller than roughly 30 grid cells in each direction, whereas longitudinal modes (i.e. shocks) may still be well resolved. As a result, the transverse kinetic energy is underestimated for up to the resolution limit . However, the plateau of almost constant for indicates that the discretisation had no significant influence on scales with wavenumbers . The effect of underestimating the transverse kinetic energy due to the discretisation of fluid variables is also observed in the ZEUS3D simulations by Pavlovski et al. (2006, Fig. 2) for wavenumbers at numerical resolution of 256^{3} grid cells. In Appendix C, we furthermore demonstrate that our results for the Fourier spectra are not affected by the specific choice of parameters of the numerical scheme for wavenumbers k 40.
6.2 Logarithmic density Fourier spectra
In analogy to the velocity Fourier spectra E(k),
we define logarithmic density fluctuation spectra
We subtract the mean logarithmic density prior to the Fourier transformation such that S(k) is a measure of density fluctuations as a function of scale. Therefore, integrating S(k) over all scales yields the square of the logarithmic density dispersion
Furthermore, integrating S(k) over the wavenumber range [k_{1},k_{2}] yields the typical density fluctuations on length scales corresponding to this range of scales.
Figure 15 (left) shows the logarithmic density fluctuation spectra S(k) together with the total velocity Fourier spectra E(k) in one plot. In contrast to the scaling of the velocity E(k), the scaling of is significantly different for solenoidal ( 0.05) and compressive forcing ( 0.09) in the inertial range.
7 variance of the velocity and density
The variance
technique provides a complementary method for measuring the scaling
exponent of Fourier spectra in the physical domain using a wavelet
transformation (Stutzki
et al. 1998). We apply the variance to our simulation
data using the tool developed and provided by Ossenkopf et al.
(2008a). This tool implements an improved version of the
original variance
(Bensch
et al. 2001; Stutzki et al. 1998).
The variance
measures the amount of structure on a given length scale
by filtering the dataset
with an
updownfunction
(typically a Frenchhat or
Mexicanhat filter) of size ,
and computing the variance of the filtered dataset. The variance is
defined as
(37) 
where the average is computed over all data points at positions . The operator stands for the convolution. We use the original Frenchhat filter with a diameter ratio of 3.0 as in previous studies using the variance technique (e.g., Ossenkopf et al. 2001; Stutzki et al. 1998; Ossenkopf & Mac Low 2002; Mac Low & Ossenkopf 2000; Ossenkopf et al. 2006).
Figure 15 (right panel) shows that the inertial range scaling obtained with the variance technique is in very good agreement with the scaling measured in the Fourier spectra. Note that the scaling exponents of Fourier spectra are ideally related to the scaling exponents of the variance by (Stutzki et al. 1998). The small deviations from this analytical relation are caused by the finite size of the dataset, the resampling procedure prior to the variance analysis applied here and the choice of the filter function (Ossenkopf et al. 2008a). However, these deviations are on the order of 4% and therefore smaller than the average snapshottosnapshot variations.
For the variance of the velocity field, , we find scaling exponents 0.05 for solenoidal forcing and 0.05 for compressive forcing. This translates into sizelinewidth relations with scaling exponents . Thus, we find 0.03 for solenoidal forcing and 0.03 for compressive forcing. Ossenkopf & Mac Low (2002) found a common powerlaw slope 0.04 for the Polaris Flare, ranging over three orders of magnitude in length scale from about down to roughly . This scaling exponent is roughly consistent with both our solenoidal and compressive forcing data, but slightly better consistent with compressive forcing. Note that the centroid velocity analysis by Ossenkopf & Mac Low (2002) is also subject to the combined effects of projectionsmoothing and intensityweighting discussed in Brunt & Mac Low (2004) and discussed in Sect. 6.1. Thus, the comparison of 3D scaling of the velocity with 2D observations should always be made with the caution that projectionsmoothing and intensityweighting roughly cancel each other out in a statistical sense (Brunt & Mac Low 2004).
We are not aware of any observational study considering the scaling of logarithmic intensity. The use of logarithmic density is useful in isothermal simulations, because the equations of hydrodynamics, Eqs. (2) and (3), are invariant under transformations in . In observations however, the intensity, T is measured instead of the density, but the intensity can be transformed into , which gives a normalised quantity similar to . This enables a straightforward comparison of simulation and observational data (yet with the limitations listed in Sect. 9). It is also interesting to look at logarithmic density and intensity scaling, because this scaling parameter is used in analytic models of the mass distribution of cores and stars by Hennebelle & Chabrier (2009,2008).
Unlike a logarithmic scaling analysis, the scaling of the linear integrated intensity, was analysed by Stutzki et al. (1998) and Bensch et al. (2001). They found for the Polaris Flare, in good agreement with the scaling exponent 0.1 obtained from solenoidal forcing in Federrath et al. (2009). In contrast, the scaling exponent obtained for compressive forcing is significantly higher ( 0.3). Bensch et al. (2001) measured scaling exponents in smallscale maps ( ) of the Polaris Flare and in Perseus/NGC 1333, which are consistent with our estimates for compressive forcing (Federrath et al. 2009, Table 1). Since both solenoidal and compressive forcings display strong intermittency at short lags (see Fig. 10), intermittency appears to be primarily measurable on scales smaller than the turbulence injection scale. Taking together the results by Bensch et al. (2001) with ours for solenoidal and compressive forcing indicates that interstellar turbulence is driven primarily on large scales, potentially with a significant amount of compressive modes present on the forcing scale (see also Brunt et al. 2009).
8 The sonic scale
The velocity Fourier spectra E(k)
discussed in Sect. 6.1 can be
described as power laws
with negative powerlaw exponents, .
This means that the typical velocity fluctuations are decreasing when
going to smaller scales. The value of the integral
over a finite range of wavenumbers with k_{1}
as the lower bound and the cutoff wavenumber
as the upper bound therefore becomes smaller with increasing k_{1}.
Thus, the turbulent flow is expected to change from a supersonic to a
subsonic flow on a certain length scale. This scale separates the
supersonic regime on large scales, where the velocity fluctuations are
supersonic from the subsonic regime, which is located on smaller
scales, where the typical velocity fluctuations are small compared to
the thermal motions of the gas. This transition scale is called the sonic
scale
.
Following Schmidt
et al. (2009), the corresponding sonic
wavenumber
in Fourier space is defined by solving the equation
implicitly for . The sonic scale is thus defined as the scale on which the mean square velocity fluctuations become comparable to the mean square of the sound speed.
We solved Eq. (38) for the sonic wavenumbers for both the solenoidal and compressive forcing cases. The sonic wavenumbers for solenoidal and compressive forcings are indicted in Fig. 15 (left) as vertical dashed lines. We find for solenoidal forcing and for compressive forcing. The corresponding sonic scales are also indicated in Fig. 15 (right) as vertical dashed lines.
The Fourier spectra S(k) shown in Fig. 15 (left) and the corresponding variance curves shown in Fig. 15 (right) for solenoidal and compressive forcing cross each other roughly at the sonic wavenumber and on the sonic scale, respectively. For compressive forcing S(k) is significantly steeper on scales larger than the sonic scale ( ) compared to scales . S(k) for compressive forcing approaches the shallower slope of S(k) for solenoidal forcing at . For there are neither significant differences between the density spectra S(k) nor the velocity spectra E(k) for solenoidal and compressive forcings.
The strong break in the logarithmic density fluctuation
spectra S(k) for
compressive forcing around
appears to be linked to the transition from supersonic motions on large
scales to subsonic motions on scales smaller than the sonic scale. In
order to quantify this, we estimated the typical density fluctuations
on supersonic scales (
)
by evaluating
.
We obtain 1.22
for solenoidal and
3.05
for compressive forcing, which is on the order of the logarithmic
density dispersions
found from the density PDFs (see Table 1). This means
that most of the power in density fluctuations is located on scales
larger than the sonic scale. In contrast, on scales smaller
than the sonic scale the typical density fluctuations can be estimated
by solving
.
We obtain 0.45
for both types of forcing. This shows that density fluctuations on
scales below the sonic scale are small compared to the typical density
fluctuations in the supersonic regime at
(see also VázquezSemadeni
et al. 2003). Moreover, Fig. 15 shows that the
typical logarithmic density fluctuations are similar for both
solenoidal and compressive forcings on scales smaller than the sonic
scale. Note that the sum of logarithmic density fluctuations on all
scales is
(39) 
for solenoidal forcing and
(40) 
for compressive forcing. As expected from Eq. (36), these values are in excellent agreement with the total logarithmic density dispersions , obtained from the density PDFs shown in Table 1.
A spatial representation of the structures exhibiting subsonic velocity dispersions is shown in Fig. 16 (bottom panel). These structures are identified in slices through the local Mach number M as regions with . Figure 16 (top panel) displays the corresponding density slices. The densityMach number correlations are quite weak, as expected for isothermal turbulence (cf. Sect. 3.6). However, Fig. 5 shows that highdensity regions exhibit lower Mach numbers on average. In real molecular clouds, the sonic scale is expected to be located on length scales within factors of a few (e.g., Goodman et al. 1998; Barranco & Goodman 1998; Falgarone et al. 1992; Schnee et al. 2007). For instance, Heyer et al. (2006) found for the Rosette MC and for G2162.5. Furthermore, the sonic scale may be associated with the transition to coherent cores (Goodman et al. 1998; Klessen et al. 2005; BallesterosParedes et al. 2003). Recent simulations of turbulent core formation by Smith et al. (2009) also suggest that starforming cores typically exhibit transonic to subsonic velocity dispersions. This can be understood if cores form close to the sonic scale in a globally supersonic turbulent medium. Figure 16 suggests that regions with subsonic velocity dispersions have different shapes and sizes for both solenoidal and compressive forcings. The movie (online version) shows that these structures are transient objects, forming and dissolving in the turbulent flow (e.g., see also VázquezSemadeni et al. 2005). If we had included selfgravity in the present study, some of these regions would have likely collapsed gravitationally, because turbulent support becomes insufficient in some of these subsonic cores (e.g., Mac Low & Klessen 2004).
Figure 16: zslices through the local density ( top panels) and Mach number fields ( bottom panels) at z=0 and t=2 T for solenoidal forcing ( left), and compressive forcing ( right). Regions with subsonic velocity dispersions ( < 1) are distinguished from regions with supersonic velocity dispersions ( > 1) in the colour scheme. The correlation between density and Mach number is quite weak. However, as shown in Fig. 5, highdensity regions exhibit lower Mach numbers on average. Thus, dense cores might naturally exhibit transonic to subsonic velocity dispersions, because their sizes are expected to be comparable to the sonic scale. The sonic scale may be the transition scale to coherent cores (e.g., Goodman et al. 1998). Although many of these ``cores'' here are transient, some of them are dense enough to become gravitationally bound, and accumulate enough mass to decouple from the overall supersonic turbulent flow. See the online version for a movie, showing the time evolution of this figure. 

Open with DEXTER 
9 Limitations
As a result of the simplicity of the hydrodynamic simulations presented in this paper, comparisons with observational data are limited and should be considered with caution. These limitations are listed below:
 We assume an isothermal equation of state, so our models are strictly speaking only applicable to molecular gas of low enough density to be optically thin to dust cooling. Variations in the equation of state can lead to changes in the density statistics (e.g., Audit & Hennebelle 2010; Li et al. 2003; Passot & VázquezSemadeni 1998). The results of the present study apply primarily to the dense interstellar molecular gas for which an isothermal equation of state is an adequate approximation (Ferrière 2001; Wolfire et al. 1995; Glover et al. 2010; Pavlovski et al. 2006).
 The numerical resolution of our simulations is limited. As shown in Fig. 6, the highdensity tails of the PDFs systematically shift to higher densities (see also Kitsionas et al. 2009; Glover et al. 2010; Price & Federrath 2010; Hennebelle & Audit 2007). However, the mean and the dispersions are well converged at the numerical resolutions of 256^{3}, 512^{3} and 1024^{3} grid points used in this study. The inertial scaling range is very small even at resolutions of 1024^{3} grid cells. However, the systematic difference in the inertial range scaling between resolutions of 512^{3} and 1024^{3} grid points is less than 3% (see Appendix C), which is less than the typical temporal variations between different realisations of the turbulent velocity and density fields.
 Our simulations adopt periodic boundary conditions. This implies that our simulations can only be representative of a subpart of a molecular cloud, for which we study turbulence statistics with highresolution numerical experiments. However, we cannot take account of the boundary effects in real molecular clouds. Simulations of largescale colliding flows (e.g., Heitsch et al. 2006; Banerjee et al. 2009; VázquezSemadeni et al. 2006; Hennebelle et al. 2008) are more suitable for studying the boundary effects during the formation of molecular clouds.
 We only analysed driven turbulence. However, there is ongoing debate about whether turbulence is driven or decaying (e.g., Stone et al. 1998; Mac Low 1999; Offner et al. 2008; Lemaster & Stone 2008). We are aware of the possibility that turbulence may in fact be excited on scales larger than the size of molecular clouds (e.g., Brunt et al. 2009), but may be globally decaying (if not replenished by a mechanism acting on galactic scales). As discussed in Sect. 2.1, this largescale decay can however act as an effective turbulence forcing on smaller scales, because kinetic energy is transported from large to small scales through the turbulence cascade.
 Centroid velocity and principal component analysis were applied to PPV cubes constructed from the simulated velocity and density fields assuming optically thin radiation transfer to estimate the intensity of emission lines. This approximation will of course not hold for optically thick tracers. A full radiative transfer calculation taking account of the level population (e.g., Steinacker et al. 2006; Keto et al. 2004; Baron et al. 2009; Pinte et al. 2009; Hauschildt & Baron 2009) of selfconsistently formed and evolved chemical tracer molecules (e.g., Glover & Mac Low 2007a; Glover et al. 2010; Glover & Mac Low 2007b) would be needed to advance on this issue.
 We neglected magnetic fields. In order to test the role of magnetic fields in star formation (e.g., Crutcher et al. 2009; Lunttila et al. 2008), we would have to include the effects of magnetic fields and ambipolar diffusion. For instance, the IMF model by Padoan & Nordlund (2002) requires magnetic fields to explain the presentday mass function, while it is still not clear whether magnetic fields are dynamically important for typical molecular clouds. However, Heyer et al. (2008) showed that magnetohydrodynamic turbulence in the Taurus MC may lead to an alignment of flows along the field lines.
 The present study did not include the effects of selfgravity, because we specifically focus on the pure turbulence statistics obtained in solenoidal and compressive forcings. In a followup study, we will include selfgravity and sink particles (e.g., Jappsen et al. 2005; Krumholz et al. 2004; Federrath et al. 2010; Bate et al. 1995) to study the influence of the different forcings on the mass distributions of sink particles. First results indicate that the sink particle formation rate is at least one order of magnitude higher for compressive forcing compared to solenoidal forcing. VázquezSemadeni et al. (2003) argue that the star formation efficiency is mainly controlled by the rms Mach number and the sonic scale of the turbulence (cf. Sect. 8). However, our preliminary results of simulations including selfgravity show that the star formation efficiency measured at a given time (i.e., the star formation rate) is much higher for compressive forcing than for solenoidal forcing with the same rms Mach number and sonic scale. This provides additional support to our main conclusion that the type of forcing must be taken into account in any theory of turbulenceregulated star formation. This needs to be investigated in future, highresolution numerical experiments including selfgravity and sink particles.
10 Summary and conclusions
We presented highresolution hydrodynamical simulations of driven isothermal supersonic turbulence, which showed that the structural characteristics of turbulence forcing significantly affect the density and velocity statistics of turbulent gas (see also Schmidt et al. 2009). We compared solenoidal (divergencefree) forcing with compressive (curlfree) turbulence forcing. Five different analysis techniques were used to compare our simulation data with existing observational data reported in the literature: probability density functions (PDFs), centroid velocity increments, principal component analysis, Fourier spectrum functions, and variances. We find that different regions in the turbulent ISM exhibit turbulence statistics consistent with different combinations of solenoidal and compressive forcing. Varying the forcing parameter in Eq. (9), we showed that a continuum of turbulence statistics exists between the two limiting cases of purely solenoidal () and purely compressive forcing (). For , turbulence behaves almost like in the case of purely solenoidal forcing, while for , turbulence is highly sensitive to changes in (cf. Fig. 8). Note that represents the natural forcing mixture used in many previous turbulence simulations. Because the behaviour of all forcing mixtures with is similar to that of purely solenoidal turbulence with (see Fig. 8), turbulence statistics is biased towards finding solenoidallike values. However, observations of regions around massive stars that drive sweptup shells into the surrounding medium (e.g., the shell in the Perseus MC and in the Rosette MC) seem better consistent with models of mainly compressive forcing (). Note that expanding HII regions around massive stars, and supernova explosions typically create such sweptup shells, which are considered to be important drivers of interstellar turbulence (Breitschwerdt et al. 2009; Mac Low & Klessen 2004)^{}. A detailed list of our results is provided below:
 1.
 The standard deviation (dispersion) of the probability distribution function (PDF) of the gas density is roughly three times larger for compressive forcing than for solenoidal forcing. This holds for both the 3D density distributions (Fig. 4 and Table 1) and the 2D column density distributions (Fig. 7 and Table 3). We extended the density dispersionMach number relations, Eqs. (18) and (19) originally investigated by Padoan et al. (1997) and Passot & VázquezSemadeni (1998). Based on the varying degree of compression obtained by solenoidal and compressive forcing, we developed a heuristic model for the proportionally constant b in the density dispersionMach number relation, which takes account of the forcing parameter (Federrath et al. 2008b). In the case of compressive forcing the proportionality constant b is close to , which confirms the result by Passot & VázquezSemadeni (1998). In contrast, solenoidal forcing yields , which is in excellent agreement with recent independent highresolution numerical simulations using solenoidal forcing (e.g., Beetz et al. 2008).
 2.
 A parameter study of eleven models with varying forcing parameter , separated by showed that the heuristic model given by Eq. (20) can only serve as a firstorder approximation to the forcing dependence of b (cf. Fig. 8). We showed that b scales with the normalised power of compressible modes in the velocity field, . A good approximation for b is given by b , where D=3 in 3D turbulence.
 3.
 We compared the density PDFs in our models with observations in the Perseus MC by Goodman et al. (2009). Goodman et al. (2009) obtained the largest density dispersion in all of the Perseus MC within a region that they call the Shell region. This Shell surrounds the massive star HD 278942, suggesting that the Shell is an expanding HII region. Sweptup shells represent geometries that can be associated with compressive turbulence forcing, because an expanding spherically symmetric shell is driven by a fully divergent velocity field. This may explain why the Shell region in the Perseus MC exhibits the largest density dispersion among all of the subregions in the Perseus MC investigated by Goodman et al. (2009). We emphasise that the Shell region does not exhibit the highest rms Mach number, but has an intermediate value among the examined subregions in the Perseus MC (Pineda et al. 2008). Furthermore, as pointed out by Goodman et al. (2009) the density dispersionMach number relation of the form given by Eq. (19) for a fixed parameter b is not observed for the Perseus MC. This apparent contradiction with Eq. (19) for a fixed parameter b is resolved, if different turbulence forcing mechanisms operate in different subregions of the Perseus MC, such that b is a function of the mixture of solenoidal and compressive modes as shown in Fig. 8.
 4.
 The turbulent density PDF is a key ingredient for the analytical models of the core mass function (CMF) and the stellar initial mass function (IMF) by Padoan & Nordlund (2002) and Hennebelle & Chabrier (2009,2008), as well as for the star formation rate models by Krumholz & McKee (2005), Krumholz et al. (2009) and Padoan & Nordlund (2009), and the star formation efficiency model by Elmegreen (2008). We showed that the dispersion of the density probability distribution is not only a function of the rms Mach number, but also depends on the nature of the turbulence forcing. All the analytical models above rely on integrals over the density PDF. Since the dispersion of the density PDF is highly sensitive to the turbulence forcing, we conclude that star formation properties derived in those analytical models are strongly affected by the assumed turbulence forcing mechanism.
 5.
 The PDFs p_{s}(s) of the logarithm of the density are roughly consistent with lognormal distributions for both solenoidal and compressive forcings. However, the distributions clearly exhibit nonGaussian higherorder moments, which are associated with intermittency. Including higherorder corrections represented by skewness and kurtosis is absolutely necessary to obtain a good analytic approximation for the PDF data, because the constraints of mass conservation (Eq. (11)) and normalisation (Eq. (12)) of the PDF must always be fulfilled. Even stronger deviations from perfect lognormal distributions are expected if the gas is nonisothermal (e.g., Li et al. 2003; Passot & VázquezSemadeni 1998; Scalo et al. 1998), magnetised (e.g., Li et al. 2008) or selfgravitating (e.g., Li et al. 2004; Federrath et al. 2008a; Kainulainen et al. 2009; Klessen 2000), which often leads to exponential wings or to powerlaw tails in the PDFs.
 6.
 NonGaussian wings of the density PDFs are a signature of intermittent fluctuations, which we further investigated using centroid velocity increments (CVIs). We find strong nonGaussian signatures for small spatial lags in the PDFs of the CVIs (Fig. 9). These PDFs exhibit values of the kurtosis significantly in excess of that expected for a Gaussian (see Fig. 10). Figure 10 can be compared with HilyBlant et al. (2008, Fig. 7), who analysed CVIs in the Taurus MC and in the Polaris Flare. The values of the kurtosis measured in the Polaris Flare are consistent with exponential values ( ) for short spatial lags, which is also compatible with the results of solenoidal forcing. In contrast, compressive forcing yields values of the kurtosis twice as large at small lags, which indicates that compressive forcing exhibits stronger intermittency. The scaling of the CVI structure functions supports the conclusion that compressive forcing exhibits stronger intermittency compared to solenoidal forcing (see Fig. 12 and Table 4). The scaling exponents of the CVI structure functions obtained for solenoidal forcing are in good agreement with the results by HilyBlant et al. (2008) obtained in the Polaris Flare for the CVI structure functions up to the 6th order using the extended selfsimilarity hypothesis.
 7.
 We applied principal component analysis (PCA) to our models. A comparison of the PCA scaling exponents with the PCA study in the Rosette MC and in G2162.5 by Heyer et al. (2006) showed that solenoidal forcing is consistent with the PCA scaling measured in G2162.5 and with the PCA scaling measured outside the HII region (Zone II) surrounding the OB star cluster NGC 2244 in the Rosette MC. On the other hand, the PCA scaling inside this HII region (Zone I) is in good agreement with the PCA scaling obtained for compressive forcing (Table 5). Similar to the Shell region in the Perseus MC, the HII region in the Rosette MC (Zone I) displays signatures of mainly compressive forcing. Recent numerical simulations by Gritschneder et al. (2009) also show that ionisation fronts driven by massive stars can efficiently excite compressible modes in the velocity field.
 8.
 The Fourier spectra of the velocity fluctuations showed that they follow power laws in the inertial range with for solenoidal forcing and for compressive forcing. Both types of forcing are therefore compatible with the scaling of velocity fluctuations inferred from observations and independent numerical simulations. The Fourier spectra of the logarithmic density fluctuations scale as for solenoidal forcing and for compressive forcing in the inertial range.
 9.
 The inertial range scaling of the velocity and logarithmic density fluctuations inferred from the Fourier spectra was confirmed using the variance technique.
 10.
 We computed the sonic scale by integrating the velocity Fourier spectra. The sonic scale separates supersonic turbulent fluctuations on large scales from subsonic turbulent fluctuations on scales smaller than the sonic scale. We found a break in the density fluctuation spectrum S(k) for compressive forcing roughly located on the sonic scale. The typical density fluctuations computed by integration of S(k) over scales larger than the sonic scale are consistent with the logarithmic density dispersions derived from the probability density functions for solenoidal and compressive forcings. On the other hand, the typical density fluctuations on scales smaller than the sonic scale are significantly smaller for both forcing types, which may reflect the transition to coherent cores (e.g., Goodman et al. 1998). Indeed, observations show that cores typically have transonic to subsonic internal velocity dispersions (e.g., Benson & Myers 1989; Foster et al. 2009; Beuther & Henning 2009; Friesen et al. 2009; Kirk et al. 2007; Lada et al. 2008; André et al. 2007; WardThompson et al. 2007). This can be understood if cores form near the sonic scale at the stagnation points of shocks in a globally supersonic turbulent ISM (cf. Sect. 8).
 11.
 We found that the correlations between the local densities and the local Mach numbers are typically quite weak (Figs. 5 and 16). However, this weak correlation shows that the local Mach number M decreases with increasing density as for solenoidal forcing and for compressive forcing for densities above the mean density. This means that dense gas tends to have smaller velocity dispersions on average, consistent with observations of dense protostellar cores.
Figure 17: Top panels: same as Fig. 15. Middle panels: same as top panels, but instead of the Fourier spectra and variances of v, the Fourier spectra and variances of the densityweighted velocity are shown. The quantity has physical reference to kinetic energy. Bottom panels: same as middle panels, but the Fourier spectra and variances of the densityweighted velocity are shown. The quantity has physical reference to a constant kinetic energy dissipation within the inertial range (Schmidt et al. 2008; Kritsuk et al. 2007). 

Open with DEXTER 
We thank Robi Banerjee, Paul Clark, Simon Glover, Alyssa Goodman, Patrick Hennebelle, Alexei Kritsuk, Alex Lazarian, Daniel Price, Stefan Schmeja, and Nicola Schneider for interesting discussions and valuable comments on the present work. We thank the referee, Chris Brunt for suggesting a parameter study with different forcing ratios , and for clarifying the effects of projectionsmoothing and intensityweighting in observations of centroid velocity maps. The variance tool used in this study was provided by Volker Ossenkopf and parallelised by Philipp Grothaus. We are grateful to Alyssa Goodman, Jaime Pineda, and Nicola Schneider for sending us their Perseus MC, and Cygnus X raw data. C. F. acknowledges financial support by the International Max Planck Research School for Astronomy and Cosmic Physics (IMPRSA) and the Heidelberg Graduate School of Fundamental Physics (HGSFP). The HGSFP is funded by the Excellence Initiative of the German Research Foundation DFG GSC 129/1. This work was partly finished while C.F. was visiting the American Museum of Natural History as a Kade fellow. R.S.K. and C.F. acknowledge financial support from the German Bundesministerium für Bildung und Forschung via the ASTRONET project STAR FORMAT (grant 05A09VHA). R.S.K. furthermore acknowledges financial support from the Deutsche Forschungsgemeinschaft (DFG) under grants No. KL 1358/1, KL 1358/4, KL 1359/5, KL 1359/10, and KL 1359/11. R.S.K. thanks for subsidies from a Frontier grant of Heidelberg University sponsored by the German Excellence Initiative and for support from the Landesstiftung BadenWürttemberg via their program International Collaboration II (grant PLSSPII/18 ). M.M. M. L. acknowledges partial support for his work from NASA Origins of Solar Systems grant NNX07AI74G. The simulations used computational resources from the HLRBII project grant h0972 at Leibniz Rechenzentrum Garching. The software used in this work was in part developed by the DOEsupported ASC/Alliance Center for Astrophysical Thermonuclear Flashes at the University of Chicago.
Appendix A: Fourier spectra and variance scaling of the combined quantities v and v
In this section we present the Fourier spectra and variance
results for the combined quantities
and
.
Usually, the pure velocity scaling is considered without density
weighting. However, for highly supersonic turbulence it is interesting
to investigate the scaling of combinations of density and velocity.
Note that CVIs (Sect. 4)
and PCA (Sect. 5)
also analyse convolutions of density and velocity statistics.
Figure A.1
(top panel) shows a repetition of Fig. 15 (scaling
of v) together with the scaling
of
(middle panel) and
(bottom panel) for direct
comparison. Since Fourier spectra and variance
analyses always represent the mean squares of these quantities,
corresponds
to the scaling of the kinetic energy density .
As shown by Kritsuk
et al. (2007) (see also Fleck 1996; Henriksen
1991), corresponds
to a constant energy flux within the inertial range. This idea
was first proposed by Lighthill
(1955). Using the eddy turnover time
as the typical evolution timescale of a turbulent fluctuation on
scale ,
the constancy of energy flux in the inertial range is
defined as
(A.1) 
which leads to the original Kolmogorov (1941) scaling (but now including density variations),
(A.2) 
for the quantity . Using the extended selfsimilarity hypothesis (Benzi et al. 1993), we showed in Schmidt et al. (2008) that the relative scaling exponents of provide a more universal scaling compared to the velocity scaling without density weighting. Figure A.1 (bottom panel) shows that the absolute scaling inferred from the Fourier spectra of is indeed close to the Kolmogorov (1941) scaling (scaling proportional to k^{5/3}) for solenoidal forcing, which is in agreement with the results obtained in Kritsuk et al. (2007). However, compressive forcing yields significantly steeper scaling (also for ), which is close to Burgers (1948) turbulence (scaling proportional to k^{2}). The corresponding results inferred from the variance analyses are compatible with the Fourier spectra to within the uncertainties. Both quantities and show breaks in the scaling close to the sonic wavenumber for compressive forcing.
Appendix B: Convergence test for the structure functions of centroid velocity increments
For an accurate and reliable determination of the structure function scaling, it must be verified that the number of data pairs used for sampling the structure functions was high enough to yield converged results. There is no general rule to determine a priori the number of data pairs necessary, because the required number of data pairs depends on the underlying statistics of the measured variable itself and on the desired structure function order. However, convergence can be tested by increasing the number of data pairs used for computing the structure functions. Showing that the computed structure functions do not change significantly by further increasing the number of data pairs demonstrates convergence. Furthermore, if convergence is verified for the highest order under consideration, then the structure functions of lower order are also converged. This is because the higherorder structure functions of a variable q reflect the statistics of higher powers of q than the lower order structure functions. This is reflected in the definition of the pth order structure function in Eq. (26).
Figure B.1 demonstrates convergence for the structure functions of CVIs with orders discussed in Sect. 4.2. We only show the compressive forcing case for clarity, but we also verified convergence for the solenoidal forcing case with the same method. Figure B.1 shows that sampling each structure function with roughly 1.7 10^{10} data pairs is sufficient to yield converged results. The total number of data pairs used to construct the CVI structure functions shown in Figs. 11 and 12 was thus roughly 81 3 1.7 10^{12} from averaging over 81 realisations of the turbulence and three projections along the x, y and zaxes for each of these realisations.
Figure B.1: The 1st (p=1) and 6th (p=6) order structure functions of the centroid velocity increments sampled with different numbers of data pairs is shown for a single snapshot at time t=2 T in zprojection for the case of compressive forcing. The number of data pairs used for sampling is given in brackets. The structure functions of centroid velocity increments are statistically converged for for sample sizes of at least 1.7 10^{10} data pairs per turbulent realisation and per projection as used throughout this study. 

Open with DEXTER 
Appendix C: Resolution study of the Fourier spectra and their dependence on the numerical scheme
The resolution and type of numerical method adopted to model supersonic turbulence are expected to critically affect the scaling of Fourier spectrum functions in the inertial range (e.g., Padoan et al. 2007; Klein et al. 2007; Kritsuk et al. 2007). In this section, we investigate the dependence of our Fourier spectra on the numerical resolution and on the numerical scheme used in the present study.
Figure C.1: Timeaveraged velocity Fourier spectra E(k) defined in Eq. (32) for numerical resolutions of 256^{3}, 512^{3} and 1024^{3} grid points obtained with solenoidal forcing ( left) and compressive forcing ( right). The inferred inertial range scaling is converged to within less than 3% at the typical resolution of 1024^{3} grid points used throughout this study for both types of forcing. 

Open with DEXTER 
Figure C.2: Dependence of the timeaveraged velocity Fourier spectra E(k) on parameters of the piecewise parabolic method (PPM) (Colella & Woodward 1984) at fixed resolution of 512^{3} grid cells. Varying the PPM diffusion parameter K between 0.0, 0.1 and 0.2 affects the dissipation range at wavenumbers . However, the effect of varying the PPM diffusion parameter is negligible for . Switching off the PPM steepening algorithm for contact discontinuities has also virtually no effect on the Fourier spectra at . 

Open with DEXTER 
C.1 Resolution study
Figure C.1 shows velocity Fourier spectra E(k) defined in Eq. (32) for numerical resolutions of 256^{3}, 512^{3} and 1024^{3} grid points. The inertial range scaling is indeed affected by the numerical resolution. For solenoidal forcing, the inertial range scaling exponent at resolution of 256^{3} grid cells is roughly 13% higher than the scaling exponent at a resolution of 1024^{3}. However, the difference between the inertial range scaling at 512^{3} and 1024^{3} is less than 3% for solenoidal forcing. For compressive forcing, the difference between the inertial range scaling exponents at resolutions of 512^{3} and 1024^{3} grid cells is less than 1%. This result indicates that the systematic dependence of the inertial range scaling on the numerical resolution is less than 3% for both solenoidal and compressive forcings. It should be emphasised that variance effects introduced by different realisations of the turbulence are typically on the order of 510% (see error bars in Fig. 15), which is higher than the systematic errors introduced by resolution effects, as long as the numerical resolution is at least 512^{3} grid cells.
C.2 Dependence on parameters of the piecewise parabolic method
We used the piecewise parabolic method (PPM) (Colella & Woodward 1984) to integrate the equations of hydrodynamics (Eqs. (2) and (3)). PPM improves on the finitevolume scheme originally developed by Godunov (1959) by representing the flow variables with piecewise parabolic functions, which makes the PPM secondorder accurate in smooth flows. However, PPM is also particularly suitable for the accurate modelling of turbulent flows involving sharp discontinuities, such as shocks and contact discontinuities. For that purpose, PPM uses a lower artificial viscosity controlled by the PPM diffusion parameter K. In three simulations with resolutions of 512^{3} grid cells, we varied the PPM diffusion parameter K between 0.0, 0.1 and 0.2. Note that K=0.1 is the value recommended by Colella & Woodward (1984), which was used for all production runs throughout this study. The PPM algorithm furthermore includes a steepening mechanism to keep contact discontinuities from spreading over too many cells. In one additional run at 512^{3}, we switched off the PPM steepening algorithm to check its influence on our results.
Figure C.2 shows that the velocity spectra E(k) decrease faster with increasing diffusion parameter K for wavenumbers . It is expected that the scheme dissipates more kinetic energy on small scales with increasing K, because the PPM diffusion algorithm is designed to act on shocks only (Colella & Woodward 1984, eq. 4.5). In contrast, Fig. C.2 demonstrates that the Fourier spectra at wavenumbers are hardly affected by the PPM diffusion algorithm for both solenoidal and compressive forcings. Note that Kritsuk et al. (2007) reported that their results for the inertial range scaling are highly sensitive to the choice of PPM diffusion parameter in the ENZO code. However, our results demonstrate that the choice of PPM diffusion parameter only affects the inertial range scaling within less than 1%, which is clearly less than the influence of the numerical resolution and less than the typical snapshottosnapshot variations. Figure C.2 furthermore demonstrates that the PPM contact discontinuity steepening has negligible effects for simulations of supersonic turbulence.
The results obtained here support our conclusion in Sect. 6 that the Fourier spectra at resolutions of 1024^{3} grid cells are robust for wavenumbers .
References
 André, P., Belloche, A., Motte, F., & Peretto, N. 2007, A&A, 472, 519
 Anselmet, F., Gagne, Y., Hopfinger, E. J., & Antonia, R. A. 1984, J. Fluid Mech., 140, 63
 Audit, E., & Hennebelle, P. 2010, A&A, 511, A76
 Azzalini, A. 1985, Scand. J. Statist., 12, 171
 BallesterosParedes, J., Klessen, R. S., & VázquezSemadeni, E. 2003, ApJ, 592, 188
 BallesterosParedes, J., Gazol, A., Kim, J., et al. 2006, ApJ, 637, 384
 Banerjee, R., VázquezSemadeni, E., Hennebelle, P., & Klessen, R. S. 2009, MNRAS, 398, 1082
 Baron, E., Hauschildt, P. H., & Chen, B. 2009, A&A, 498, 987
 Barranco, J. A., & Goodman, A. A. 1998, ApJ, 504, 207
 Bate, M. R., Bonnell, I. A., & Price, N. M. 1995, MNRAS, 277, 362
 Beetz, C., Schwarz, C., Dreher, J., & Grauer, R. 2008, Phys. Lett. A, 372, 3037
 Bensch, F., Stutzki, J., & Ossenkopf, V. 2001, A&A, 366, 636
 Benson, P. J., & Myers, P. C. 1989, ApJS, 71, 89
 Benzi, R., Ciliberto, S., Tripiccione, R., et al. 1993, Phys. Rev. E, 48, 29
 Beuther, H., & Henning, T. 2009, A&A, 503, 859
 Boldyrev, S. 2002, ApJ, 569, 841
 Boldyrev, S., Nordlund, Å., & Padoan, P. 2002, ApJ, 573, 678
 Breitschwerdt, D., de Avillez, M. A., Fuchs, B., & Dettbarn, C. 2009, Space Sci. Rev., 143, 263
 Brunt, C. M., & Heyer, M. H. 2002a, ApJ, 566, 276
 Brunt, C. M., & Heyer, M. H. 2002b, ApJ, 566, 289
 Brunt, C. M., & Mac Low, M. 2004, ApJ, 604, 196
 Brunt, C. M., Heyer, M. H., VázquezSemadeni, E., & Pichardo, B. 2003, ApJ, 595, 824
 Brunt, C. M., Heyer, M. H., & Mac Low, M. 2009, A&A, 504, 883
 Brunt, C. M., Federrath, C., & Price, D. J. 2010, MNRAS, 403, 1507
 Burgers, J. M. 1948, Adv. Appl. Mech., 1, 171
 Burkhart, B., FalcetaGonçalves, D., Kowal, G., & Lazarian, A. 2009, ApJ, 693, 250
 Colella, P., & Woodward, P. R. 1984, J. Comput. Phys., 54, 174
 Crutcher, R. M., Hakobian, N., & Troland, T. H. 2009, ApJ, 692, 844
 Dib, S., Brandenburg, A., Kim, J., Gopinathan, M., & André, P. 2008, ApJ, 678, L105
 Dobler, W., Haugen, N. E. L., Yousef, T. A., & Brandenburg, A. 2003, Phys. Rev. E, 68, 026304
 Dubey, A., Fisher, R., Graziani, C., et al. 2008, in Numerical Modeling of Space Plasma Flows, ed. N. V. Pogorelov, E. Audit, & G. P. Zank, ASP Conf. Ser., 385, 145
 Dubinski, J., Narayan, R., & Phillips, T. G. 1995, ApJ, 448, 226
 Elmegreen, B. G. 2002, ApJ, 577, 206
 Elmegreen, B. G. 2008, ApJ, 672, 1006
 Elmegreen, B. G., & Lada, C. J. 1977, ApJ, 214, 725
 Elmegreen, B. G., & Scalo, J. 2004, ARA&A, 42, 211
 Esquivel, A., Lazarian, A., Horibe, S., et al. 2007, MNRAS, 381, 1733
 Eswaran, V., & Pope, S. B. 1988, Comput. Fluids, 16, 257
 Falgarone, E., & Phillips, T. G. 1990, ApJ, 359, 344
 Falgarone, E., Puget, J.L., & Perault, M. 1992, A&A, 257, 715
 Falgarone, E., Lis, D. C., Phillips, T. G., et al. 1994, ApJ, 436, 728
 Federrath, C., Glover, S. C. O., Klessen, R. S., & Schmidt, W. 2008a, Phys. Scr. T, 132, 014025
 Federrath, C., Klessen, R. S., & Schmidt, W. 2008b, ApJ, 688, L79
 Federrath, C., Klessen, R. S., & Schmidt, W. 2009, ApJ, 692, 364
 Federrath, C., Banerjee, R., Clark, P. C., & Klessen, R. S. 2010, ApJ, 713, 269
 Ferrière, K. M. 2001, Rev. Mod. Phys., 73, 1031
 Fleck, Jr., R. C. 1996, ApJ, 458, 739
 Foster, J. B., Rosolowsky, E. W., Kauffmann, J., et al. 2009, ApJ, 696, 298
 Friesen, R. K., Di Francesco, J., Shirley, Y. L., & Myers, P. C. 2009, ApJ, 697, 1457
 Frisch, U. 1995, Turbulence (Cambridge Univ. Press)
 Fryxell, B., Olson, K., Ricker, P., et al. 2000, ApJS, 131, 273
 Glover, S. C. O., & Mac Low, M.M. 2007a, ApJS, 169, 239
 Glover, S. C. O., & Mac Low, M.M. 2007b, ApJ, 659, 1317
 Glover, S. C. O., Federrath, C., Mac Low, M.M., & Klessen, R. S. 2010, MNRAS, in press [arXiv:0907.4081]
 Godunov, S. K. 1959, Math. Sbornik, 47, 271
 Goodman, A. A., Barranco, J. A., Wilner, D. J., & Heyer, M. H. 1998, ApJ, 504, 223
 Goodman, A. A., Pineda, J. E., & Schnee, S. L. 2009, ApJ, 692, 91
 Gritschneder, M., Naab, T., Walch, S., Burkert, A., & Heitsch, F. 2009, ApJ, 694, L26
 Haugen, N. E., & Brandenburg, A. 2004, Phys. Rev. E, 70, 026405
 Hauschildt, P. H., & Baron, E. 2009, A&A, 498, 981
 Heitsch, F., Mac Low, M.M., & Klessen, R. S. 2001, ApJ, 547, 280
 Heitsch, F., Slyz, A. D., Devriendt, J. E. G., Hartmann, L. W., & Burkert, A. 2006, ApJ, 648, 1052
 Hennebelle, P., & Audit, E. 2007, A&A, 465, 431
 Hennebelle, P., & Chabrier, G. 2008, ApJ, 684, 395
 Hennebelle, P. & Chabrier, G. 2009, ApJ, 702, 1428
 Hennebelle, P., Banerjee, R., VázquezSemadeni, E., Klessen, R. S., & Audit, E. 2008, A&A, 486, L43
 Henriksen, R. N. 1991, ApJ, 377, 500
 Heyer, M. H., & Brunt, C. M. 2004, ApJ, 615, L45
 Heyer, M. H., & Schloerb, F. P. 1997, ApJ, 475, 173
 Heyer, M. H., Williams, J. P., & Brunt, C. M. 2006, ApJ, 643, 956
 Heyer, M., Gong, H., Ostriker, E., & Brunt, C. 2008, ApJ, 680, 420
 Heyer, M., Krawczyk, C., Duval, J., & Jackson, J. M. 2009, ApJ, 699, 1092
 HilyBlant, P., Falgarone, E., & Pety, J. 2008, A&A, 481, 367
 Jappsen, A.K., Klessen, R. S., Larson, R. B., Li, Y., & Mac Low, M.M. 2005, A&A, 435, 611
 Kainulainen, J., Beuther, H., Henning, T., & Plume, R. 2009, A&A, 508, L35
 Keto, E., Rybicki, G. B., Bergin, E. A., & Plume, R. 2004, ApJ, 613, 355
 Kirk, H., Johnstone, D., & Tafalla, M. 2007, ApJ, 668, 1042
 Kissmann, R., Kleimann, J., Fichtner, H., & Grauer, R. 2008, MNRAS, 391, 1577
 Kitsionas, S., Federrath, C., Klessen, R. S., et al. 2009, A&A, 508, 541
 Klein, R. I., Inutsuka, S.I., Padoan, P., & Tomisaka, K. 2007, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil, 99
 Klessen, R. S. 2000, ApJ, 535, 869
 Klessen, R. S. 2001, ApJ, 556, 837
 Klessen, R. S., Heitsch, F., & Mac Low, M.M. 2000, ApJ, 535, 887
 Klessen, R. S., BallesterosParedes, J., VázquezSemadeni, E., & DuránRojas, C. 2005, ApJ, 620, 786
 Kolmogorov, A. N. 1941, Dokl. Akad. Nauk SSSR, 32, 16
 Kowal, G., Lazarian, A., & Beresnyak, A. 2007, ApJ, 658, 423
 Kravtsov, A. V. 2003, ApJ, 590, L1
 Kritsuk, A. G., Norman, M. L., Padoan, P., & Wagner, R. 2007, ApJ, 665, 416
 Krumholz, M. R., & McKee, C. F. 2005, ApJ, 630, 250
 Krumholz, M. R., McKee, C. F., & Klein, R. I. 2004, ApJ, 611, 399
 Krumholz, M. R., McKee, C. F., & Tumlinson, J. 2009, ApJ, 699, 850
 Lada, C. J., Muench, A. A., Rathborne, J., Alves, J. F., & Lombardi, M. 2008, ApJ, 672, 410
 Larson, R. B. 1969, MNRAS, 145, 271
 Larson, R. B. 1981, MNRAS, 194, 809
 Larson, R. B. 2005, MNRAS, 359, 211
 Lazarian, A., & Esquivel, A. 2003, ApJ, 592, L37
 Lemaster, M. N., & Stone, J. M. 2008, ApJ, 682, L97
 Lemaster, M. N., & Stone, J. M. 2009, ApJ, 691, 1092
 Li, Y., Klessen, R. S., & Mac Low, M.M. 2003, ApJ, 592, 975
 Li, P. S., Norman, M. L., Mac Low, M.M., & Heitsch, F. 2004, ApJ, 605, 800
 Li, P. S., McKee, C. F., Klein, R. I., & Fisher, R. T. 2008, ApJ, 684, 380
 Lighthill, M. J. 1955, in Gas Dynamics of Cosmic Clouds, IAU Symp., 2, 121
 Lis, D. C., Pety, J., Phillips, T. G., & Falgarone, E. 1996, ApJ, 463, 623
 Lis, D. C., Keene, J., Li, Y., Phillips, T. G., & Pety, J. 1998, ApJ, 504, 889
 Lombardi, M., Alves, J., & Lada, C. J. 2006, A&A, 454, 781
 Lunttila, T., Padoan, P., Juvela, M., & Nordlund, Å. 2008, ApJ, 686, L91
 Mac Low, M.M. 1999, ApJ, 524, 169
 Mac Low, M.M., & Klessen, R. S. 2004, Rev. Mod. Phys., 76, 125
 Mac Low, M.M., & Ossenkopf, V. 2000, A&A, 353, 339
 Mac Low, M.M., Klessen, R. S., Burkert, A., & Smith, M. D. 1998, Phys. Rev. Lett., 80, 2754
 Mac Low, M.M., Balsara, D. S., Kim, J., & de Avillez, M. A. 2005, ApJ, 626, 864
 McKee, C. F., & Ostriker, E. C. 2007, ARA&A, 45, 565
 Miesch, M. S., & Bally, J. 1994, ApJ, 429, 645
 Miesch, M. S., Scalo, J., & Bally, J. 1999, ApJ, 524, 895
 Moisy, F., & Jiménez, J. 2004, J. Fluid Mech., 513, 111
 Murtagh, F., & Heck, A. 1987, Multivariate Data Analysis, Astrophys. Space Sci. Library, 131
 Myers, P. C. 1983, ApJ, 270, 105
 Nordlund, Å., & Padoan, P. 1999, in Interstellar Turbulence, ed. J. Franco, & A. Carraminana, 218
 Offner, S. S. R., Klein, R. I., & McKee, C. F. 2008, ApJ, 686, 1174
 Ossenkopf, V., & Mac Low, M.M. 2002, A&A, 390, 307
 Ossenkopf, V., Klessen, R. S., & Heitsch, F. 2001, A&A, 379, 1005
 Ossenkopf, V., Esquivel, A., Lazarian, A., & Stutzki, J. 2006, A&A, 452, 223
 Ossenkopf, V., Krips, M., & Stutzki, J. 2008a, A&A, 485, 917
 Ossenkopf, V., Krips, M., & Stutzki, J. 2008b, A&A, 485, 719
 Ostriker, E. C., Gammie, C. F., & Stone, J. M. 1999, ApJ, 513, 259
 Ostriker, E. C., Stone, J. M., & Gammie, C. F. 2001, ApJ, 546, 980
 Padoan, P., & Nordlund, Å. 2002, ApJ, 576, 870
 Padoan, P., & Nordlund, A. 2009, ApJ, submitted [arXiv:0907.0248]
 Padoan, P., Nordlund, Å., & Jones, B. J. T. 1997, MNRAS, 288, 145
 Padoan, P., Boldyrev, S., Langer, W., & Nordlund, Å. 2003, ApJ, 583, 308
 Padoan, P., Jimenez, R., Nordlund, Å., & Boldyrev, S. 2004, Phys. Rev. Lett., 92, 191102
 Padoan, P., Juvela, M., Kritsuk, A., & Norman, M. L. 2006, ApJ, 653, L125
 Padoan, P., Nordlund, Å., Kritsuk, A. G., Norman, M. L., & Li, P. S. 2007, ApJ, 661, 972
 Passot, T., & VázquezSemadeni, E. 1998, Phys. Rev. E, 58, 4501
 Pavlovski, G., Smith, M. D., & Mac Low, M.M. 2006, MNRAS, 368, 943
 Penston, M. V. 1969, MNRAS, 144, 425
 Perault, M., Falgarone, E., & Puget, J. L. 1986, A&A, 157, 139
 Pineda, J. E., Caselli, P., & Goodman, A. A. 2008, ApJ, 679, 481
 Pinte, C., Harries, T. J., Min, M., et al. 2009, A&A, 498, 967
 Porter, D. H., Pouquet, A., & Woodward, P. R. 1992a, Theor. Comput. Fluid Dynam., 4, 13
 Porter, D. H., Pouquet, A., & Woodward, P. R. 1992b, Phys. Rev. Lett., 68, 3156
 Porter, D. H., Pouquet, A., & Woodward, P. R. 1994, Phys. Fluids, 6, 2133
 Press, W. H., Flannery, B. P., Teukolsky, S. A., & Vetterling, W. T. 1986, Numerical recipes (Cambridge University Press), 818 S.
 Price, D. J., & Federrath, C. 2010, MNRAS, submitted
 Ridge, N. A., Schnee, S. L., Goodman, A. A., & Foster, J. B. 2006, ApJ, 643, 932
 Sánchez, N., Alfaro, E. J., & Pérez, E. 2005, ApJ, 625, 849
 Scalo, J., & Elmegreen, B. G. 2004, ARA&A, 42, 275
 Scalo, J., VázquezSemadeni, E., Chappell, D., & Passot, T. 1998, ApJ, 504, 835
 Schmidt, W., Hillebrandt, W., & Niemeyer, J. C. 2006, Comput. Fluids, 35, 353
 Schmidt, W., Federrath, C., & Klessen, R. 2008, Phys. Rev. Lett., 101, 194505
 Schmidt, W., Federrath, C., Hupp, M., Kern, S., & Niemeyer, J. C. 2009, A&A, 494, 127
 Schnee, S., Caselli, P., Goodman, A., et al. 2007, ApJ, 671, 1839
 She, Z.S., & Leveque, E. 1994, Phys. Rev. Lett., 72, 336
 Smith, R. J., Clark, P. C., & Bonnell, I. A. 2009, MNRAS, 396, 830
 Solomon, P. M., Rivolo, A. R., Barrett, J., & Yahil, A. 1987, ApJ, 319, 730
 Steinacker, J., Bacmann, A., & Henning, T. 2006, ApJ, 645, 920
 Stone, J. M., Ostriker, E. C., & Gammie, C. F. 1998, ApJ, 508, L99
 Stutzki, J., Bensch, F., Heithausen, A., Ossenkopf, V., & Zielinsky, M. 1998, A&A, 336, 697
 Tamburro, D., Rix, H.W., Leroy, A. K., et al. 2009, AJ, 137, 4424
 Tassis, K. 2007, MNRAS, 382, 1317
 VázquezSemadeni, E. 1994, ApJ, 423, 681
 VázquezSemadeni, E., BallesterosParedes, J., & Klessen, R. S. 2003, ApJ, 585, L131
 VázquezSemadeni, E., Kim, J., Shadmehri, M., & BallesterosParedes, J. 2005, ApJ, 618, 344
 VázquezSemadeni, E., Ryu, D., Passot, T., González, R. F., & Gazol, A. 2006, ApJ, 643, 245
 Vincent, A., & Meneguzzi, M. 1991, J. Fluid Mech., 225, 1
 Wang, J., Townsley, L. K., Feigelson, E. D., et al. 2008, ApJ, 675, 464
 Wang, J., Feigelson, E. D., Townsley, L. K., et al. 2009, ApJ, 696, 47
 WardThompson, D., André, P., Crutcher, R., et al. 2007, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil, 33
 Wolfire, M. G., Hollenbach, D., McKee, C. F., Tielens, A. G. G. M., & Bakes, E. L. O. 1995, ApJ, 443, 152
Online Material
Sol_vs_Comp_slices_FederrathEtAl2010.wmv
Footnotes
 ... forcing^{}
 A movie is only available in electronic form at http://www.aanda.org
 ... ^{}
 Note that Passot & VázquezSemadeni (1998) suffers from a number of typographical errors as a result of lastminute change of notation. Please see Mac Low et al. (2005, footnote 5) for a number of corrections.
 ... range^{}
 By its formal definition for incompressible turbulence studies (e.g., Frisch 1995), the inertial range is the range of scales for which the turbulence statistics are not directly influenced by the forcing acting on scales larger than the inertial range, and not directly influenced by the viscosity acting on scales smaller than the inertial range. The inertial range is typically very small in numerical experiments, because of the high numerical viscosity caused by the discretisation scheme, given the resolutions achievable with current computer technology (see also Sect. 6 and Appendix C).
 ... NGC 2244^{}
 The formation of the star cluster XA in the Rosette MC was likely triggered by the accumulation of material in the expanding shell surrounding the OB star cluster NGC 2244 (Wang et al. 2009,2008). This emphasises the importance of expanding HII regions in triggering subsequent star formation by compression of gas in expanding shells (Elmegreen & Lada 1977).
 ... ^{}
 See also Tamburro et al. (2009) for an observational study.
All Tables
Table 1: Statistical moments and fit parameters of the PDFs of the volumetric density for solenoidal and compressive forcing shown in Fig. 4.
Table 2: Standard deviations of the density PDFs as a function of numerical resolution for solenoidal and compressive forcing shown in Fig. 6.
Table 3: Same as Table 1, but for the PDFs of the column density shown in Fig. 7.
Table 4: Scaling of the structure functions of centroid velocity increments.
Table 5: Comparison of measured PCA scaling slopes.
All Figures
Figure 1: Ratio of compressive power to the total power in the turbulence force field. The solid lines labelled with 1D, 2D, and 3D show the analytical expectation for this ratio, Eq. (9), as a function of the forcing parameter for one, two and threedimensional forcing, respectively. The diamonds and squares show results of numerical simulations in 3D and 2D with , separated by . Those models were run at a numerical resolution of 256^{3} and 1024^{2} grid points in 3D and 2D, respectively. The two extreme forcing cases of purely solenoidal forcing () and purely compressive forcing () are indicated as ``sol'' and ``comp'', respectively. Note that in any 1D model, all power is in the compressive component, and thus , independent of . 

Open with DEXTER  
In the text 
Figure 2: Maps showing density ( top), vorticity ( middle) and divergence ( bottom) in projection along the zaxis at time t=2 T as an example for the regime of statistically fully developed, compressible turbulence for solenoidal forcing ( left) and compressive forcing ( right). Top panels: column density fields in units of the mean column density. Both maps show three orders of magnitude in column density with the same scaling and magnitudes for direct comparison. Middle panels: projections of the modulus of the vorticity . Regions of intense vorticity appear to be elongated filamentary structures often coinciding with positions of intersecting shocks. Bottom panels: projections of the divergence of the velocity field showing the positions of shocks. Negative divergence corresponds to compression, while positive divergence corresponds to rarefaction. 

Open with DEXTER  
In the text 
Figure 3: Top panel: minimum and maximum logarithmic density as a function of the dynamical time T. Note that compressive forcing yields much stronger compression and rarefaction compared to solenoidal forcing, although the rms Mach number is roughly the same in both cases (see Federrath et al. 2009, Fig. 2). Bottom panel: rms vorticity and rms divergence as a function of the dynamical time. Within the first 2 T, a statistically steady state was reached for both solenoidal (sol) and compressive (comp) forcing. This allows us to average statistical measures (probability density functions, centroid velocity increments, principal component analysis, Fourier spectra and variances) in the range to improve statistical significance of our results and to estimate the amplitude of temporal fluctuations (snapshottosnapshot variations) between different realisations of the turbulence. 

Open with DEXTER  
In the text 
Figure 4: Volumeweighted density PDFs p(s) of the logarithmic density in linear scaling ( top panel), which displays the peak best, and in logarithmic scaling ( bottom panel) to depict the low and highdensity wings. The PDF obtained from compressive forcing ( ) is significantly wider than the solenoidal one ( ). The peak is shifted to lower values of the logarithmic density s, because of mass conservation, defined in Eq. (11). The density PDF from solenoidal forcing is compatible with a Gaussian distribution. However, there are also nonGaussian features present, which are associated with intermittency effects. These are more prominent in the density PDF obtained from compressive forcing, exhibiting statistically significant deviations from a perfect lognormal (fit using Eq. (10) shown as dashed lines). A skewed lognormal fit (dashdotted lines) given by Eq. (14) provides a better representation, but still does not fit the highdensity tail of the PDF obtained for compressive forcing. Both the PDF data obtained from solenoidal and compressive forcing are best described as lognormal distributions with higherorder corrections defined in Eq. (17), which take into account both the nonGaussian skewness and kurtosis of the distributions. These fits are shown as solid lines (skewkurtlognormal fit). The first four standardised moments defined in Eqs. (13) of the distributions in and s are summarised in Table 1 together with the fit parameters. The grey shaded regions indicate 1 error bars due to temporal fluctuations of the distributions in the regime of fully developed, supersonic turbulence. A total number of 1024^{3} data points contribute to each PDF. 

Open with DEXTER  
In the text 
Figure 5: Volumeweighted correlation PDFs of local Mach number M versus logarithmic density s for solenoidal ( left) and compressive forcing ( right). Adjacent contour levels are spaced by in probability density. Density and Mach number exhibit a very weak, but nonzero correlation in both forcing cases, which provides an explanation for the nonGaussian features in the density PDFs of Fig. 4 (VázquezSemadeni 1994; Passot & VázquezSemadeni 1998; Kritsuk et al. 2007). The two solid lines, intersecting the maxima of both distributions show the mean Mach number as a function of the logarithmic density . The tendency for highdensity gas having lower Mach numbers on average is indicated as power laws in the highdensity parts of the distributions. This suggests that the Mach number for solenoidal and for compressive forcing. 

Open with DEXTER  
In the text 
Figure 6: Density PDFs at numerical resolutions of 256^{3}, 512^{3} and 1024^{3} grid cells. The PDFs show very good overall convergence, especially around the peaks. Table 2 shows that the standard deviations are converged with numerical resolution. The highdensity tails, however, are not converged even at a numerical resolution of 1024^{3} grid points, indicating a systematic shift to higher densities with resolution. This limitation is shared among all turbulence simulations (see also, Kitsionas et al. 2009; Price & Federrath 2010; Hennebelle & Audit 2007). The lowdensity wings are subject to strong temporal fluctuations due to intermittent bursts caused by headon collisions of shocks followed by strong rarefaction waves (e.g., Kritsuk et al. 2007). The intermittency causes deviations from a perfect Gaussian distribution and accounts for nonGaussian higherorder moments (skewness and kurtosis) in the distributions. 

Open with DEXTER  
In the text 
Figure 7: Same as Fig. 4, but the time and projectionaveraged logarithmic column density PDFs of are shown. and denote the column density and the mean column density integrated along the i=x, y, z principal axes respectively. As for the volumetric PDFs of Fig. 4, the standard deviation of the column density PDF obtained from compressive forcing is roughly three times larger than from solenoidal forcing (see Table 3). The inset in the upper right corner shows the PDFs of column density computed in zprojection at a fixed time t=2 T, corresponding to the snapshots shown in Fig. 2. The density dispersions computed for these instantaneous PDFs are and for solenoidal forcing, and and for compressive forcing. Although these distributions are quite noisy, the influence of the forcing is still clearly discernible. Thus, by studying instantaneous column density PDFs, which are accessible to observations, one should be able to distinguish solenoidal from compressive forcing. 

Open with DEXTER  
In the text 
Figure 8: Diamonds: the proportionality parameter b in the density dispersionMach number relation, Eq. (18), computed as for eleven 3D models at numerical resolution of 256^{3} grid cells ( top panel) and eleven 2D models at numerical resolution of 1024^{2} grid cells ( bottom panel), ranging from purely compressive forcing () to purely solenoidal forcing (). The parameter b decreases smoothly from for compressive forcing to in 3D and in 2D for solenoidal forcing. Stars: ratio of longitudinal to total power in the velocity power spectrum (see Sect. 6.1). This quantity provides a measure for the relative amount of compression induced by the turbulent velocity field, and appears to be correlated with the standard deviation of the density PDF. Squares: same as stars, but multiplied by the geometrical factor with D=3 for the threedimensional case and D=2 for the twodimensional case. The quantity provides a good numerical estimate of the PDF parameter b. The dashed lines show model fits using Eq. (23) for D=3 ( top panel) and D=2 ( bottom panel). 

Open with DEXTER  
In the text 
Figure 9: PDFs of centroid velocity increments, computed using Eqs. (24) and (25) are shown as a function of lag in units of grid cells for solenoidal forcing ( left) and compressive forcing ( right). The PDFs are very close to Gaussian distributions for long lags, whereas for short lags, they develop exponential tails, which is a manifestation of intermittency (e.g., HilyBlant et al. 2008, and references therein). 

Open with DEXTER  
In the text 
Figure 10: Kurtosis of the PDFs of centroid velocity increments shown in Fig. 9 as a function of the lag in units of grid cells for solenoidal and compressive forcing. Note that a kurtosis value of 3 (horizontal dotdashed line) corresponds to the value for a Gaussian distribution. NonGaussian values of the kurtosis are obtained for . The error bars contain both snapshottosnapshot variations as well as the variations between centroid velocity increments computed by integration along the x, y and z axes. This figure can be compared to observations of the Polaris Flare and Taurus MC (see Fig. 7 of HilyBlant et al. 2008). 

Open with DEXTER  
In the text 
Figure 11: Scaling of the structure functions of centroid velocity increments defined in Eq. (26) for solenoidal forcing ( left) and compressive forcing ( right) up to the 6th order. Scaling exponents obtained using powerlaw fits following Eq. (27) within the inertial range are indicated in the figures and summarised in Table 4. 

Open with DEXTER  
In the text 
Figure 12: Same as Fig. 11, but using the extended selfsimilarity hypothesis (Benzi et al. 1993), allowing for a direct comparison of the scaling exponents of centroid velocity increments with the study by HilyBlant et al. (2008) for the Polaris Flare and Taurus MC (see Table 4). 

Open with DEXTER  
In the text 
Figure 13: Principal component analysis (PCA) for solenoidal ( left) and compressive forcing ( right). The PCA slopes obtained for solenoidal and compressive forcings are summarised and compared with observations by Heyer et al. (2006) in Table 5. The error bars contain the contribution from temporal variations and from three different projections along the x, y and zaxes. The data were resampled from 1024^{3} to 256^{3} grid points prior to PCA. The resampling speeds up the PCA and has virtually no effect on the inertial range scaling (see e.g., Federrath et al. 2009; Padoan et al. 2006). 

Open with DEXTER  
In the text 
Figure 14: Top panels: total, transverse (rotational) and longitudinal (compressible) velocity Fourier spectra E(k) defined in Eq. (32) and compensated by k^{2} for solenoidal ( left) and compressive forcing ( right). Error bars indicate temporal variations, which account for an uncertainty of roughly of all scaling slopes reported for the inertial range . The inferred inertial range scaling exponents for both solenoidal and compressive forcing are consistent with independent numerical simulations and with observations of the sizelinewidth relation (see text). Note that the transverse part, falls off more steeply than the longitudinal part, for both forcing types in the inertial range. Bottom panels: ratio of the energy in longitudinal velocity modes to the total energy in velocity modes . For solenoidal forcing, we obtain 1/3 in the inertial range (horizontal dashdotted line), because compression can only occur in one of the three spatial dimensions on average (Federrath et al. 2008b; Elmegreen & Scalo 2004). For compressive forcing, this ratio is roughly 1/2, which corresponds to an equipartition of longitudinal and transverse velocity modes. Note however that compressive forcing can compress the gas in all three spatial dimensions directly, whereas solenoidal forcing can only induce compression indirectly through the velocity field (Federrath et al. 2008b). The excess of longitudinal modes at high wavenumbers stems from numerical dissipation, which is more effectively dissipating transverse than longitudinal modes on small scales due to the discretisation onto a grid. This suggests that roughly 30 grid cells are needed to accurately resolve a vortex, while a shock is typically resolved with roughly 3 grid cells using the piecewise parabolic method (Colella & Woodward 1984). However, for a numerical resolution of 1024^{3} grid cells, we find that wavenumbers are almost unaffected by the discretisation and by the parameters of the numerical scheme (see Appendix C). 

Open with DEXTER  
In the text 
Figure 15: Left panel: fourier spectra of the velocity, E(k) defined in Eq. (32) (crosses and diamonds) and Fourier spectra of the logarithmic density fluctuations, S(k) defined in Eq. (35) (triangles and squares) for solenoidal and compressive forcing, respectively. Both E(k) and S(k) are compensated by k^{2} allowing for a better determination of the inertial range scaling. The density fluctuation power spectra differ significantly in the inertial range with for solenoidal and for compressive forcing. The scale on which the density fluctuation spectra from solenoidal and compressive forcing cross each other and where the slope obtained in compressive forcing breaks and approaches the shallower slope of the solenoidal forcing case roughly coincides with the sonic wavenumber (vertical dashed lines) defined in Eq. (38). Right panel: same as left panel, but instead of using Fourier spectra to determine the inertial range scaling, we use the variance method to derive the scaling slopes in physical space. Note that the scaling slopes obtained with the variance technique are related to the slopes of the Fourier spectra by (Stutzki et al. 1998). Error bars denote 1 temporal fluctuations. 

Open with DEXTER  
In the text 
Figure 16: zslices through the local density ( top panels) and Mach number fields ( bottom panels) at z=0 and t=2 T for solenoidal forcing ( left), and compressive forcing ( right). Regions with subsonic velocity dispersions ( < 1) are distinguished from regions with supersonic velocity dispersions ( > 1) in the colour scheme. The correlation between density and Mach number is quite weak. However, as shown in Fig. 5, highdensity regions exhibit lower Mach numbers on average. Thus, dense cores might naturally exhibit transonic to subsonic velocity dispersions, because their sizes are expected to be comparable to the sonic scale. The sonic scale may be the transition scale to coherent cores (e.g., Goodman et al. 1998). Although many of these ``cores'' here are transient, some of them are dense enough to become gravitationally bound, and accumulate enough mass to decouple from the overall supersonic turbulent flow. See the online version for a movie, showing the time evolution of this figure. 

Open with DEXTER  
In the text 
Figure 17: Top panels: same as Fig. 15. Middle panels: same as top panels, but instead of the Fourier spectra and variances of v, the Fourier spectra and variances of the densityweighted velocity are shown. The quantity has physical reference to kinetic energy. Bottom panels: same as middle panels, but the Fourier spectra and variances of the densityweighted velocity are shown. The quantity has physical reference to a constant kinetic energy dissipation within the inertial range (Schmidt et al. 2008; Kritsuk et al. 2007). 

Open with DEXTER  
In the text 
Figure B.1: The 1st (p=1) and 6th (p=6) order structure functions of the centroid velocity increments sampled with different numbers of data pairs is shown for a single snapshot at time t=2 T in zprojection for the case of compressive forcing. The number of data pairs used for sampling is given in brackets. The structure functions of centroid velocity increments are statistically converged for for sample sizes of at least 1.7 10^{10} data pairs per turbulent realisation and per projection as used throughout this study. 

Open with DEXTER  
In the text 
Figure C.1: Timeaveraged velocity Fourier spectra E(k) defined in Eq. (32) for numerical resolutions of 256^{3}, 512^{3} and 1024^{3} grid points obtained with solenoidal forcing ( left) and compressive forcing ( right). The inferred inertial range scaling is converged to within less than 3% at the typical resolution of 1024^{3} grid points used throughout this study for both types of forcing. 

Open with DEXTER  
In the text 
Figure C.2: Dependence of the timeaveraged velocity Fourier spectra E(k) on parameters of the piecewise parabolic method (PPM) (Colella & Woodward 1984) at fixed resolution of 512^{3} grid cells. Varying the PPM diffusion parameter K between 0.0, 0.1 and 0.2 affects the dissipation range at wavenumbers . However, the effect of varying the PPM diffusion parameter is negligible for . Switching off the PPM steepening algorithm for contact discontinuities has also virtually no effect on the Fourier spectra at . 

Open with DEXTER  
In the text 
Copyright ESO 2010