Issue 
A&A
Volume 630, October 2019



Article Number  A97  
Number of page(s)  21  
Section  Interstellar and circumstellar matter  
DOI  https://doi.org/10.1051/00046361/201833970  
Published online  26 September 2019 
How do velocity structure functions trace gas dynamics in simulated molecular clouds?
^{1}
MaxPlanckInstitut für Astronomie,
Königstuhl 17,
69117
Heidelberg,
Germany
email: roxanaadela.chira@alumni.uniheidelberg.de
^{2}
I. Physikalisches Institut, Universität zu Köln,
Zülpicher Straße 77,
50937
Köln,
Germany
email: ibanez@ph1.unikoeln.de
^{3}
MaxPlanckInstitut für Extraterrestrische Physik,
Giessenbachstrasse 1,
85748
Garching,
Germany
^{4}
Department of Astrophysics, American Museum of Natural History,
79th St. at Central Park West,
New York,
NY
10024,
USA
email: mordecai@amnh.org
^{5}
Zentrum für Astronomie, Institut für Theoretische Astrophysik, Universität Heidelberg,
AlbertUeberleStr. 2,
69120
Heidelberg,
Germany
^{6}
Center for Computational Astrophysics, Flatiron Institute,
162 Fifth Ave,
New York,
NY
10010,
USA
Received:
27
July
2018
Accepted:
9
August
2019
Context. Supersonic disordered flows accompany the formation and evolution of molecular clouds (MCs). It has been argued that this is turbulence that can support against gravitational collapse and form hierarchical substructures.
Aims. We examine the time evolution of simulated MCs to investigate: What physical process dominates the driving of turbulent flows? How can these flows be characterised? Are they consistent with uniform turbulence or gravitational collapse? Do the simulated flows agree with observations?
Methods. We analysed three MCs that have formed selfconsistently within kiloparsecscale numerical simulations of the interstellar medium (ISM). The simulated ISM evolves under the influence of physical processes including selfgravity, stratification, magnetic fields, supernovadriven turbulence, and radiative heating and cooling. We characterise the flows using velocity structure functions (VSFs) with and without density weighting or a density cutoff, and computed in one or three dimensions. However, we do not include optical depth effects that can hide motions in the densest gas, limiting comparison of our results with observations.
Results. In regions with sufficient resolution, the densityweighted VSFs initially appear to follow the expectations for uniform turbulence, with a firstorder powerlaw exponent consistent with Larson’s sizevelocity relationship. Supernova blast wave impacts on MCs produce shortlived coherent motions at large scales, increasing the scaling exponents for a crossing time. Gravitational contraction drives smallscale motions, producing scaling coefficients that drop or even turn negative as small scales become dominant. Removing the density weighting eliminates this effect as it emphasises the diffuse ISM.
Conclusions. We conclude that two different effects coincidentally reproduce Larson’s size velocity relationship. Initially, uniform turbulence dominates, so the energy cascade produces VSFs that are consistent with Larson’s relationship. Later, contraction dominates and the densityweighted VSFs become much shallower or even inverted, but the relationship of the global average velocity dispersion of the MCs to their radius follows Larson’s relationship, reflecting virial equilibrium or freefall collapse. The injection of energy by shocks is visible in the VSFs, but decays within a crossing time.
Key words: turbulence / ISM: kinematics and dynamics / ISM: structure / ISM: clouds
© R.A. Chira et al. 2019
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Open Access funding provided by Max Planck Society.
1 Introduction
It has long been known that star formation preferentially occurs within molecular clouds (MCs). However, the physics of the star formation process is still not completely understood. It is obvious that gravity is the key factor for star formation as it drives collapse motions and operates on all scales. However, one needs additional processes that stabilise the gas or terminate star formation quickly in order to explain the low star formation efficiencies observed in MCs. Although there are many processes that act at the different scales of MCs, turbulent support has often been argued to be the best candidate for this task.
In the literature, turbulence plays an ambiguous role in the context of star formation. In most of the cases, turbulence is expected to stabilise MCs on large scales (Fleck 1980; McKee & Zweibel 1992; Mac Low 2003), while feedback processes and shear motions heavily destabilise or even disrupt cloudlike structures (Tan et al. 2013; Miyamoto et al. 2014). However, it remains unclear whether there are particular mechanisms that dominate the driving of turbulence within MCs, as every process is supposed to be traced by typical features in the observables. Yet, these features are either not seen or are too ambiguous to clearly reflect the dominant driving mode. For example, turbulence that is driven by largescale velocity dispersions during global collapse (BallesterosParedes et al. 2011a,b; Hartmann et al. 2012) produces PCygni line profiles that have not yet been observed on scales of entire MCs. Internal feedback, on the contrary, seems more promising as it drives turbulence from small to large scales (Dekel & Krumholz 2013; Krumholz et al. 2014). Observations, though, demonstrate that the required driving sources need to act on scales of entire clouds, which typical feedback, such as radiation, winds, jets, or supernovae (SNe), cannot achieve (Heyer & Brunt 2004; Brunt et al. 2009; Brunt & Heyer 2013).
There have been many theoretical studies that have examined the nature and origin of turbulence within the various phases of the interstellar medium (ISM; Mac Low & Klessen 2004, and references within). The most established characterisation of turbulence in general was introduced by Kolmogorov (1941) who investigated fully developed, incompressible turbulence driven on scales larger than the object of interest, and dissipating on scales much smaller than those of interest. In the scope of this paper this object is a single MC. MCs are highly compressible, though. Only a few analytical studies have treated this case. She & Lévêque (1994) and Boldyrev (2002), for example, generalised and extended the predicted scaling of the decay of turbulence to supersonic turbulence. Galtier & Banerjee (2011) and Banerjee & Galtier (2013) provided an analytic description of the scaling of massweighted structure functions.
Larson (1981) found a relation between the linewidth σ and the effective radius R of MCs. Subsequent investigators have settled on the form of the relation being (Solomon et al. 1987; Falgarone et al. 2009; Heyer et al. 2009) (1)
Goodman et al. (1998) showed that analysis techniques used to study this relation could be distinguished by whether they studied single or multiple clouds using single or multiple tracer species. Explanations for this relation have relied on either turbulent cascades (Larson 1981; Kritsuk et al. 2013a, 2015; Gnedin 2015; Padoan et al. 2016), or the action of selfgravity (Elmegreen 1993, 2007; VázquezSemadeni et al. 2006; Heyer et al. 2009; BallesterosParedes et al. 2011b).
These can potentially be distinguished by examining the velocity structure function (VSF). Kritsuk et al. (2013a) carefully reviews the argument for Larson’s sizevelocity relation depending on the turbulent cascade. In short, in an energy cascade typical for turbulence, the secondorder structure function has a lag dependence ℓ^{ζ(2)} with ζ(2) ≃ 1∕2. In IbáñezMejía et al. (2016, hereafter Paper I) the authors argue that uniform driven turbulence was unable to explain the observed relation in a heterogeneous ISM, but that the relation could be naturally explained by hierarchical gravitational collapse.
In this paper, we examine the velocity structure functions of three MCs that formed selfconsistently from SNdriven turbulence in the simulations by Paper I and IbáñezMejía et al. (2017, hereafter Paper II). We study how the turbulence within the clouds’ gas evolves. The key questions we address are the following: What dominates the turbulence within the simulated MCs? Does the observed linewidthsize relation arise from the turbulent flow? How can structure functions inform us about the evolutionary state of MCs and the relative importance of largescale turbulence, discrete blast waves, and gravitational collapse?
In Sect. 2, we introduce the simulated clouds in the context of the underlying physics involved in the simulations. Furthermore, we describe the theoretical basics of velocity structure functions. Section 3 demonstrates that the velocity structure function is a useful tool to characterise the dominant driving mechanisms of turbulence in MCs and can be applied to both simulated and observed data. We examine the influence of using onedimensional velocity measurements, different Jeans refinement levels, density thresholds, density weighting on the applicability of the velocity structure function, and the results obtained with it in Sect. 4. At the end of that section, we also compare our results to observational studies. We summarise our findings and conclusions in Sect. 5. The simulation data and the scripts that this work is based on are published in the Digital Repository of the American Museum of Natural History (Chira et al. 2018a).
2 Methods
2.1 Cloud models
The analysis in this paper is based on a sample of three MCs identified within a threedimensional (3D), magnetohydrodynamical, adaptive mesh refinement simulation using the FLASH code (Fryxell et al. 2000). Paper I and Paper II, as well as Chira et al. (2018b, Paper III hereafter), describe the simulations and the clouds in more detail. We summarise the most relevant properties.
The numerical simulation models a 1 × 1 × 40 kpc^{3} section of the multiphase, turbulent ISM of a disc galaxy, where dense structures form selfconsistently in convergent, turbulent flows Paper I. The model includes gravity – a background galacticdisc potential accounting for a stellar component and a dark matter halo, as well as selfgravity turned on after 250 Myr of simulated time – SNdriven turbulence, photoelectric heating and radiative cooling, and magnetic fields. Although hundreds of dense clouds form within the simulated volume, Paper II focussed on three clouds, which were resimulated at a much higher spatial resolution. The internal structures of the MCs are resolved using adaptive mesh refinement, focussing grid resolution on dense regions where Jeans unstable structures must be resolved with a minimum of 4 cells (λ_{J} > 4Δx_{min}). For a maximum resolution of Δx = 0.1 pc, the corresponding maximum resolved density is 8 × 10^{3} cm^{−3} for gas at a temperature of 10 K (e.g. Paper II, Eq. (15)). We define MCs as regions above a fixed number density threshold with fiducial value n_{cloud} = 100 cm^{−3}. We chose this threshold as it roughly corresponds to the density when CO becomes detectable. The MCs have initial masses at the onset of selfgravity of 3 × 10^{3}, 4 × 10^{3}, and 8 × 10^{3} M_{⊙} and are denoted as M3, M4, and M8, respectively, hereafter. In this paper, we use the data within (40 pc)^{3} subregions centred on the highresolution clouds’ centres of mass, which we map to a uniform grid of 0.1 pc zones for analysis. For illustrations of the morphologies of the three clouds we refer to Fig. 1 of Paper III.
It is important to point out that the clouds are embedded within a complex turbulent environment, gaining and losing mass as they evolve. Paper II described the time evolution of the properties of all three clouds in detail, in particular, mass, size, velocity dispersion, and accretion rates, in the context of MC formation and evolution within a galactic environment. Paper III studied the properties, evolution, and fragmentation of filaments that selfconsistently condense within these clouds. We paid particular attention to the accuracy of typical stability criteria for filaments, comparing the results to the theoretical predictions, showing that simplified analytic models do not capture the complexity of fragmentation due to their simplifying assumptions.
2.2 Velocity structure function
In this paper, we probe the power distribution of turbulence throughout the entire simulated MCs by using the velocity structure function (VSF). The VSF is a twopoint correlation function, (2)
that measures the mean velocity difference (3)
between two points x and x + ℓ, with ℓ being the direction vector pointing from the first to the second point. The VSF S_{p} is usually reported as a function of lag distance, ℓ = ℓ, between the correlated points. The coherent velocity differences measured by the VSF can be produced by both the energy cascade expected in turbulent flows, and by coherent motions such as collapse, rotation, or blast waves. Those patterns become more prominent the higher the value of the power p is (Heyer & Brunt 2004).
For fully developed, homogeneous, isotropic, turbulence the VSF is welldescribed by a powerlaw relation (Kolmogorov 1941; She & Lévêque 1994; Boldyrev 2002): (4)
Kolmogorov (1941) predicts that the thirdorder exponent, ζ(3), is equal to unity for an incompressible flow. As a consequence the kinetic energy decays with , with being the wavenumber of the turbulence mode. For a supersonic flow, however, ζ(3) > 1 is expected. Based on Kolmogorov’s work, She & Lévêque (1994) and Boldyrev (2002) extended and generalised the analysis and predicted the following intermittency corrections to Kolmogorov’s scaling law. For incompressible turbulence with filamentary dissipative structures She & Lévêque (1994) predict that the VSFs scale with power law index (5)
while supersonic flows with sheetlike dissipative structures are predicted to scale with (Boldyrev 2002) (6)
It should be noted that both equations return a value of ζ(3) = 1, but this is only an exact result for the She & Lévêque model, while it is a result of normalisation in the case of Boldyrev.
In the case of compressible turbulence, the energy cascade can no longer be expressed in terms of a pure velocity difference because density fluctuations become important. Turbulence should then show a cascade in some densityweighted VSF analogous to the incompressible case. Padoan et al. (2016) defined a densityweighted VSF to attempt to capture this process, which we use in our subsequent analysis (7)
Alternatives have been proposed by Kritsuk et al. (2013b) based on an analysis of the equations of compressible flow that should be explored in future work.
In many cases a threedimensional computation of the VSF cannot be performed because of the observational constraint that only lineofsight velocities are available. We therefore compare our threedimensional (3D) results to onedimensional (1D), densityweighted VSFs (8)
with e_{i} representing the unit vector along the i = x, y, or zaxis.
Benzi et al. (1993) introduced the principle of extended selfsimilarity. It proposes that there is a constant relationship between the scaling exponents of different orders at all lag scales so that ζ can be measured from S_{p}∕S_{3}, which typically gives a clearer powerlaw behaviour. The selfsimilarity parameter is defined as (9)
As mentioned before, both Eqs. (5) and (6) return values of ζ(3) = 1. Therefore, those equations also offer predictions for Z(p).
For the discussion below, we measure ζ(p) by fitting a powerlaw, given by (10)
with A being the proportionality factor of the powerlaw to the simulated measurements. We choose the smallest lag of the fitting range to be equal to eight zones, sufficiently large to ensure that our fit does not include the numerical dissipation range. For more details of the fitting procedure we refer to Appendix A.
We follow observational practice and reduce the computational effort of this study by generally focussing on clouds defined by a density threshold. However, Paper II shows that there is usually no sharp increase in density between the ISM and the clouds. Instead, the gas becomes continuously denser towards the centres of mass within the clouds. Consequently, our use of a density threshold is a somewhat artificial boundary between the clouds and the ISM. Observationally, however, introducing a column density (or intensity) threshold is unavoidable, be it due to technical limitations (e.g. detector sensitivity) or the nature of the underlying physical processes (for example, excitation rates, or critical densities). Therefore, we also study how a density threshold influences the VSF and its evolution.
At our fiducial density threshold, we actually consider only ≤1.5% of the volume in the high resolution cube. To understand the influence of this limitation we set up a test scenario (see Sect. 3.4) by removing the density threshold (setting n_{cloud} = 0 cm^{−3}) that results in analysing the entire data cube. Details of the method for computing the VSFs in these two cases are given in Appendix A.
As in the case without a density threshold it would be too computationally expensive to compute all lags to all zones. Thus, we randomly choose a set of 5% of the total number of zones as reference points and only compute relative velocities from the entire cube to these zones. By choosing the starting points randomly we ensure that all parts of the cubes are considered. As a consequence, there is only a small likelihood (5% × 1.5% = 0.075%) that any given zone chosen will be within the cloud. Therefore, we emphasise that it is likely that the two subsamples (no density threshold and cloudonly) do not have a common subset of starting vectors. Nevertheless, the random sample still includes > 4 × 10^{3} zones in the cloud, so the sample does include information on VSFs of material in the cloud.
3 Results
3.1 Examples
In this section, we present our results on how VSFs reflect the velocity structure within and around MCs.
Figure 1 shows nine examples of densityweighted VSFs (Eq. (7)). The figure shows the VSFs of all three clouds (columns) at times of 1.0, 3.0, and 4.2 Myr after the onset of gravity. All plots show orders p = 1–3. The solid lines show the fitted powerlaw relations as given by Eq. (10).
The examples demonstrate that, in general, the measured VSFs cannot be described by a single powerlaw relation over the entire range of ℓ. Instead they are composed of roughly three different regimes: one at small scales at 0.8 pc ≲ ℓ ≲ 3 pc, a second one within 3 pc ≲ ℓ ≲ 10–15 pc, and the last one at large scales with 10–15 pc ≲ ℓ ≲ 30 pc. We find that only the small and intermediate ranges may be represented by a common powerlaw relation. On larger scales, one observes a local minimum before the VSFs either increase or stagnate. Additional examples of VSFs are given in Appendix B.
The examples in Fig. 1 and Appendix B illustrate how VSFs react to different scenarios that affect the turbulent structure of the entire clouds. All clouds at t = 1.0 Myr show the case where turbulence is driven on large scales and naturally decays towards smaller scales. This is the most common behaviour seen in all three MCs within the first ~1.5 Myr of the simulations. During this interval of time the clouds experience the effect of selfgravity for the first time in their evolution and need to adjust to this new condition. Until this occurs, their VSFs continue to be dominated by the freely cascading turbulence that previously dominated the kinetic structure of the clouds. We note that with each refinement it takes finite time for the turbulence to propagate to smaller scale, so the cloud evolution at high resolution was started well before selfgravity was turned on (see Paper II, Seifried et al. 2017).
The later examples represent the clouds when the VSFs are dominated by sources that drive the flow within the clouds in a more coherent way. M3 at t = 3.0 Myr and M4 at t = 4.2 Myr show VSFs at times when the clouds have just been hit by a SN blast wave. One clearly sees how the amplitude of the VSFs are increased byan order of magnitude or more compared to the time before. Especially the power at small scales below a few parsecs is highly amplified as a result of the shock. Despite the increase of turbulent power at small scales, a large amount of energy is injected at large scales, as well. However, the effect of SN shocks last for only a short period of time (see Sect. 3.2).
M8 at t = 3.0 Myr demonstrates the imprint of gravitational contraction. Here, the VSF is almost flat, or even slightly increasing towards smaller separation scales. This kind of profile is typical for gas that is gravitationally contracting (Boneberg et al. 2015; Burkhart et al. 2015). Gas moves into the inner regions of the cloud, reducing the average lag distances, while being accelerated by the infall to higher velocity. As a consequence, large amounts of kinetic energy are transferred to smaller scales and higher densities, flattening the corresponding densityweighted VSF.
Fig. 1
Examples of VSFs from models (left to right) M3, M4, and M8 as function of lag scale ℓ and order p, based on data with density threshold. The examples are given for times (top to bottom) t = 1.0 Myr, 3.0 Myr, and 4.2 Myr. The dots (connected by dashed lines) show the values computed from the simulations. The solid lines represent the powerlaw relation fitted to the respective structure functions. 
3.2 Time evolution
Figure 2a summarises the time evolution of the powerlaw index ζ(p) fit to the densityweighted VSF obtained for each cloud, and each order p. The figure shows several interesting features. First, initially, at t = 0 Myr, all calculated values of ζ are above the predicted values (see Eqs. (5) and (6)). This probably occurs primarily because the base simulation only resolves down to 0.49 pc before additional grid levels are added to resolve the clouds, so it cannot resolve the turbulence inertial range below approximately 3 pc. This can be seen in the t = 0 Myr power spectrum in Fig. 25, Appendix B, of Paper II. As the zoomin simulations evolve, the turbulence cascades to smaller scales in dense regions that are better resolved, so the densityweighted VSF slopes initially drop to the values expected for supersonic turbulence. The VSFs without densityweighting and, even more so, without density threshold remain too steep despite the increased resolution in dense regions, because these VSFs remain dominated by numerical diffusion in the diffuse gas that has not been further refined.
Second, ζ for all orders decreases with time as the clouds are first refined and then begin to gravitationally collapse. Distributed gravitational collapse causes an increase in relative velocities at increasingly small scales as material falls into filaments and nodes. The increase in smallscale power leads to a flattening or even inversion of the VSF and thus a decrease in ζ. Third, occasionally one observes bumps and dips in slope in all orders of VSFs (e.g. M3 or M8 around t = 1.7 Myr). These features only last for short periods of time (up to 0.6 Myr), but set in rather abruptly and represent sudden changes in largescale power that change the VSF slope.
Figure 3a shows the corresponding time evolution of the selfsimilarity parameter, Z. One sees that most of the time the measured values of Z are in agreement or at least closely approaching the predicted values. The occasional peaks in Z (for example, in M4 at t = 4.1 Myr) occur at times when the scaling exponents of the VSFs ζ(3) reach values close to or below 0. A decrease in Z (for example, in M3 around t = 1.8 Myr) occurs whenSN shocks hit and heavily impact the clouds, producing stronger effects in higher order VSFs.
We note that the nearby SNe that we mark in Figs. 2 and 3 and consider in the analysis were listed by Paper II as exploding at a radius of up to 100 pc from the clouds’ centres of mass. The shock fronts move at average speeds of 50–100 km s^{−1} through the ISM, so it can take the blast waves more than 1 Myr to reach the clouds. Thus, it is important to keep in mind that the MCs do not react immediately to SNe, and that the time between the explosion of a SN and the interaction of its shock front with one of the clouds varies depending on the distance, as well as with the composition of the ISM along the propagation path.
Paper II discusses the influence of these SNe on the modelled clouds. In particular, the authors focus on the overall properties of the clouds, like total mass, accretion rate, total velocity dispersion, and evolution. They find that there is a tight correlation between those properties and that the clouds evolution strongly depends on whether the clouds are hit by blast waves, as well as the details of these interactions (meaning, for example, the strengthor direction of the shock jump). Therefore, it is intuitive that blast waves have a substantial influence on the turbulence within the clouds, as well. We discuss this in more detail in the following sections of this paper.
In the rest of this section, we study how VSFs computed in different ways compare to these density weighted results. We compare thefindings with the results we have obtained with our original setup. In Sect. 4, we discuss and interpret these resultsin more detail.
3.3 Lineofsight VSF
Previously, we have seen how the VSF behaves and evolves within the clouds. To do this, we derived the relative velocities based on the 3D velocity vectors from the simulations. However, observed VSFs can only be measured using lineofsight velocities. In this subsection we investigate how VSFs derived from 1D relative velocities compare to the 3D VSFs presented before.
Figures 2b and 3b show measured ζ and Z, respectively, derived based on Eq. (8). We see that, in most of the cases, most of the 1D VSFs agree well with each other, as well as with the corresponding 3D VSFs. The VSFs of M4 are the best examples for this.
However, there are also cases in which the 1D VSF temporarily evolves completely differently than the 3D VSF. For example, the 1D VSFalong the xaxis in M3 initially behaves like the corresponding 3D VSF, though with lower absolute values of ζ (or higher values of Z). However, during the period t = 2.5–3.8 Myr the functions diverge. While the 3D ζ decreases further and switches sign, the ζ based on the1D VSF along the xaxis shows a local maximum before converging with the 3D ζ again. Another example for temporal divergence are the VSFs of M8 projectedonto the x and zaxes. Here, the VSFs develop below or above the level of the 3D VSF, respectively. After the impact of the SN blast wave at t = 1.5 Myr all three VSFs converge with each other, as well as with the 3D VSF under the influence of the ongoing gravitational collapse.
Fig. 2
Time evolution of scaling exponent ζ(p) of the pth order VSF. Panels a–d: the measurements for M3 (left), M4 (middle), and M8 (right), respectively. Of these, panel a represents the standard analysis while the other rows illustrate the results of different variations in the analysis as noted in the figures. Panel e: the values of ζ measured within M3 as a function of the Jeans refinement level the cloud has been modelled with. We note that these more expensive runs were not run for as long as the fiducial run. In all panels, the grey dotted vertical lines mark the times than a SN explodes within 100 pc of the cloud, while the blue areas indicate periods of enhanced mass accretion onto the clouds. The coloured horizontal lines show the predicted values for incompressible turbulence (dashdotted lines; She & Lévêque 1994) and for highly compressible, supersonic turbulence (dashed lines; Boldyrev 2002). 
3.4 Density thresholds
We now examine the VSFs of the entire data cubes without setting a density threshold (i.e. setting n_{cloud} = 0 cm^{−3}). Figure 2c shows ζ, while Fig. 3c shows Z in this case.These figures clearly illustrate that the measurements in the samples without density threshold completely differ from those with the density threshold. The measured values of ζ are much higher in the ISM than in the cloudonly sample. Although we see a slight decline of ζ in M4 and M8 as the gas contracts under the influence of gravity in the vicinity of the clouds, ζ generally evolves differently here than in the analysis that is focussed on the clouds, remaining very steep. We see a high rate of random fluctuations in the evolution of ζ, as well. Furthermore, contrary to all of our other test scenarios, we see here that all Z are constant in time and within all clouds, with values slightly lower than those predicted by She & Lévêque (1994) for incompressible flows. This is consistent with the high sound speed in the hot gas that fills most of the volume of the computational box, which results in subsonic flows predominating.
Fig. 3
Like Fig. 2, but for the measured selfsimilarity parameter Z = ζ(p)∕ζ(3) of the pth order VSF. 
3.5 Density weighting
As mentioned in Sect. 2.2, Eq. (2) represents the definition of the densityweighted VSF. While this represents the observational situation better, the theoretical predictions were developed for the unweighted statistic. Thus a comparison of results for the two variations in our model is of interest. There are a few studies that have targeted this question (e.g. Benzi et al. 1993, 2010; Schmidt et al. 2008; Gotoh et al. 2002). However, all of them considered isotropic, homeogeneous, turbulent flows that are not comparable to our clouds. Padoan et al. (2016) use both methods, but not on the same set of data.
In this section, we investigate the influence of density weighting on VSFs by repeating the original analysis with the nonweighted VSF given by Eq. (2). Figures 2d and 3d show the measured values of ζ and Z derived from the nonweighted VSFs, respectively.
Comparing the weighted and nonweighted samples, we see the following: The nonweighted ζ (Fig. 2d) traces the interactions between the gas of the clouds and the SN shocks in the same way as occurs for the densityweighted VSF. In M3 and M8 we also see that the values of ζ decrease as the clouds evolve, yet not as fast as they do in the densityweighted VSFs, which focus attention on the dense regions of strongest collapse. The measurements in M4, however, are almost constant over time. In all the cases, the values of ζ never decrease below 0.5; a behaviour that clearly differs from what we have observed in the densityweighted VSFs. The density weighting weakens the influence of the highly compressible flows in the densest regions, but not so much as in the case with no density threshold. Consequently, the evolution of Zs (see Fig. 3d) becomes smoother, as well, as there is no sign inversion of ζ. As a result the values of Z fluctuate slightly when shocks hit, and otherwise vary between the compressible and incompressible limits.
3.6 Jeans length refinement
The results we have discussed so far are based on simulation data presented in Paper I and Paper II. Due to the huge computational expense of the variety of physical and numerical processes (fluid dynamics, adaptive mesh refinement, SNe, magnetic fields, radiative heating and cooling, and many more) within those simulations, though, they have required some compromises. One of these compromises was the Jeans refinement criterion used as part of the adaptive mesh refinement mechanism. The authors resolved local Jeans lengths by only four zones (λ_{J} = 4Δx). This is the minimal requirement for modelling selfgravitating gas in discs in order to avoid artificial fragmentation (Truelove et al. 1998). Other studies, have shown that a significantly higher refinement is needed to reliably resolve turbulent structures and flows within discs to resolve turbulence (Federrath et al. 2011; Turk et al. 2012). In our case, the key question is how quickly the turbulent cascade fills in after the multiple steps of refinement to higher resolution required to develop the high resolution cubes we use. Although we have a different physical situation, the earlier results still emphasise the importance of how well the Jeans length is resolved.
In the appendix of Paper II, the authors examine the effect the number of zones used for the Jeans refinement has on the measured kinetic energy. For this, they have rerun the simulations of M3 twice; once with λ_{J} = 8Δx for the first 3 Myr after selfgravity was activated, and once with λ_{J} = 32Δx for the firstmegayear of the cloud’s evolution. The authors show that the λ_{J} = 32Δx simulations smoothly recover the energy power spectrum on all scales already after this first megayear. The other two setups do this, as well, but only over longer timescales (see also Paper II, Seifried et al. 2017).
Furthermore, Paper II calculated the difference in the cloud’s total kinetic energy as a function of time and refinement level. They found that the λ_{J} = 4Δx simulations miss a significant amount of kinetic energy, namely up to 13% compared to λ_{J} = 8Δx and 33% compared to λ_{J} = 32Δx. However, they also observed that these differences peak around t = 0.5 Myr and decrease afterwards, as the λ_{J} = 4Δx and λ_{J} = 8Δx simulations adjust to the new refinement levels. Thus, the results we have derived from the λ_{J} = 4Δx simulations need to be evaluated with respect to this lack of turbulent energy, although the clouds’ dynamics remains dominated by gravitational collapse. It also means that the λ_{J} = 4Δx data become more reliable the longer the simulations evolve.
In order to study how the level of Jeans refinement influences the behaviour of the VSFs, we investigate the M3 data of the λ_{J} = 8Δx and λ_{J} = 32Δx simulations. Figures 2e and 3e show ζ and Z for λ_{J} = 8Δx and λ_{J} = 32Δx. In Fig. 4 we directly compare the measurements of all refinement levels relative to λ_{J} = 4Δx.
The λ_{J} = 8Δx model shows the same behaviour as λ_{J} = 4Δx, with values in both samples being in good agreement as the top panel of Fig. 4 demonstrates. Over the entire observed time span, the measured values of ζ decrease as the VSF become flatter. At the time the SNe interact with the cloud, over the course of about a megayear after traveling across the distance from the point of explosion to the cloud, the VSFs steeply increase toward larger scales, causing values of ζ (Fig. 2e) to jump. Compared to the λ_{J} = 4Δx sample, the peak in ζ is smoother and lasts longer at higher Jeans resolution.
These same effects can be seen in Fig. 3e where the drop of Z due to the SN shock lasts longer than it did for λ_{J} = 4Δx. Besides this, the time evolution of Z for λ_{J} = 8Δx is as sensitive to the turbulencerelated events as it was for λ_{J} = 4Δx. The divergence produced when gravity has transferred the majority of power to smaller scales occurs at the same time. The actual depth of the drop is a numerical artefact caused by ζ(3) being equal or close to zero at this very time step.
The picture changes when we analyse the VSFs based on the λ_{J} = 32Δx runs (Figs. 2e, 3e, and 4). Here one sees that the measured values of both ζ (Fig. 2e) and Z (Fig. 3e) are similar to those for λ_{J} = 4Δx for the first 0.2 Myr. After this short period, though, the evolution of ζ diverges. While ζ(1) and ζ(2) continue to decrease similar to λ_{J} = 4Δx but at lower rate, ζ(3) increases until it peaks at t = 0.8 Myr and falls steeply again. This divergence has notable impact on the evolution of Z, as well. The bottom panel of Fig. 4 illustrates the different evolution of measured ζ and Z in the two simulations more clearly. One sees that the differences between the samples follow the same pattern for all orders of p. The differences, though, increase with the order: While the values for ζ(1) are still in good agreement, the measured values of ζ(2) and ζ(3) for λ_{J} = 32Δx are 40% and 100% higher than those measured for λ_{J} = 4Δx, respectively. Consequently, this causes differences in Z(p) of 30–52% between the simulations. At t = 1.2 Myr, the last time step of this sample, the values of all ζ equal the measurements of λ_{J} = 4Δx again. As the cost of extending the λ_{J} = 32Δx simulation is prohibitive, we cannot determine whether this agreement will continue.
4 Discussion
4.1 Time evolution
We have seen in Sect. 3 that densityweighted VSFs reflect a combination of uniform, compressible turbulence, largescale shocks, and gravitational collapse. Extended selfsimilarity emphasises the turbulent nature of these highReynolds numbers flows even in regions of gravitational collapse. The measured values of ζ differ from the predicted values by She & Lévêque (1994) and Boldyrev (2002) for most of the time of the clouds’ evolution. The origin of these differences are to be found in a combination of numerical dissipation and the underlying physical conditions. Contrary to the conditions in the reference models, we examine a multiphase medium with magnetised turbulence that is hardly uniform, isothermal, or homogeneous. However, the diffuse medium is also underresolved compared to the dense gas, thanks to our Jeans refinement strategy, which removes smallscale power in the diffuse medium.
The impact of SN shocks hitting the clouds is to inject power at all scales (Fig. 1). The resulting VSFs tend to lose their powerlaw character. Fitting a powerlaw to them anyway results in substantial perturbations from the predictions for compressible turbulence even under extended selfsimilarity. Figure 3 shows times of SN explosions and periods of strong accretion onto the clouds. Remembering that it can take the shock front more than 1 Myr to propagate from the site of the SN explosion to the molecular cloud, perturbations in Z not associated with zerocrossings by ζ(3) are consistent with being caused by SN shock front interactions with the clouds. These shock interactions last for only a fraction of a megayear, though, consistent with the crossing time of the blast wave through the dense interior of the cloud, after which the turbulent nature of the flow reasserts itself.
As the clouds gravitationally collapse, the resulting increase in smallscale power flattens or even inverts the densityweighted VSFs, resulting in decreasing or even negative values of ζ (Fig. 2a). The increase in smallscale power can also be derived from the increasingly negative binding energy of the clouds as further gas falls into them Paper II. At the same time as the turbulence becomes increasingly nonuniform and anisotropic because of the importance of gravitation, the bulk velocity dispersion of the cloud increases. Paper I showed that Eq. (1) is satisfied at these late times, but not at early times, less than a freefall time, when the velocity dispersion inherited from the background turbulent flow is independent of the size of the cloud. These early times are when the turbulence dominates the flow and the secondorder power law is roughly ζ(2) ≃ 1∕2. This suggests that the apparent agreement with Larson’s sizevelocity relationship is coincidental. Observing twopoint correlations using a method that cannot capture the dense flows adequately will yield this result from the cloud envelopes, though, thus perhaps explaining the apparent success of such efforts (Goodman et al. 1998, Fig. 9 shows how these different interpretations can arise).
Extended selfsimilarity shows VSF ratios characteristic of compressible turbulence (Fig. 3a), as can be seen from their tending to lie between the incompressible limit of She & Lévêque (1994) and the extremely compressible Burgers turbulence limit of Boldyrev (2002). (The extended selfsimilarity procedure fails as ζ(3) passes through zero, however, so it must be interpreted in concert with the raw values of ζ.) This suggests that, just as the extended selfsimilarity procedure removes the effects of dissipation, it also removes the effects of hierarchical gravitational collapse, while continuing to reflect the turbulence in these high Reynoldsnumbers flows.
Fig. 4
Comparison of VSF scaling exponents, ζ (left), and selfsimilarity parameters, Z (right), depending on the Jeans refinement of the simulation runs the data are based on. The abscissas give values from λ_{J} = 4Δx, while the ordinates give values from λ_{J} = 8Δx (top) and λ_{J} = 32Δx (bottom). All data points refer to the M3 cloud and represent different lags in the same time step in the respective simulations. 
4.2 Lineofsight velocities
In Sect. 3.3 we have seen that the ζ and Z derived from the 1D VSFs generally evolve similarly to those derived from the 3D VSFs. Yet, we have also seen that individual sight lines may evolve differently. These differences appear to reflect the detailed geometry of shock impacts on the cloud, which are reflected more strongly in the higherorder VSFs. For example, for the first 2 Myr of the evolution of M4 the values of Z along the yaxis are significantly higher than those observed along the other axes and diverge significantly from the values expected for uniform turbulence. It should be recalled here a perturbation in Z usually corresponds to an episode of strong shock driving, suggesting an impact along the yaxis at this time. Along the other two axes, Z continues to agree with supersonic turbulence (Boldyrev 2002). This effect is only visible as we analyse the three dimensions separately, while the driving of the gas along the yaxis is averaged out in the 3D VSFs (see Fig. 3a).
In summary, for a fully developed 3D turbulent field we expect that 1D VSFs behave similarly to 3D VSFs. However, when there is a preferred direction along which the gas flows, the 1D and 3D VSFs differ significantly from each other. Thus, we predict that observed VSFs reflect the nature of turbulence within MCs unless there is clear evidence that the gas is driven in a particular direction (e.g. by a stellar wind or SN shock front).
We note that this analysis does not take typical lineofsight effects, such as optical depth or blending, into account. Future studies needto investigate this point in more detail by performing VSF analyses based on full radiative transfer calculations.
4.3 Density thresholds
We find that the structure and behaviour of VSFs strongly depends on whether or not we assume a density threshold in computing them. In the fiducial case, where n_{cloud} = 100 cm^{−3}, we have seena mostly straight decline of ζ while Z remains fairly constant over time, reflecting the contraction of the clouds due to selfgravity. On the other hand, if we remove the density threshold, including the entire highresolution cube in the calculation, we observe a completely different picture. The high velocities present in the diffuse ISM surrounding the cloud lead to strong largescale power and thus much steeper VSFs, corresponding to higher values of ζ. There is still a slightly declining trend in ζ, but the evolution is dominated by random fluctuations. We also see that Z remains constant for the entire simulation in every case. However, the VSF scaling exponents are about four times steeper than the values that are predicted by (Boldyrev 2002) for incompressible flows. This is partly an effect of numerical dissipation reducing the velocities in lessrefined diffuse gas. However, it is probably also due to the sharp reduction in the velocity dispersion of dense gas already noted in Paper I, given that the dense gas has small characteristic scales. This suggests that they are dominated by the subsonic flow in the hot gas with T > 10^{6} K that occupies most of the volume of the box.
Furthermore, the results demonstrate that the effect of SNe, and the interaction of the produced shock fronts with the ISM acts rather locally. This means that a single SN can not drive the gas dynamics on scales of our entire cubes (40^{3} pc^{3}), at least not in the same way as it does on scales of individual MCs. Rather the VSFs reflect the superposition of multiple SNe that only combined drive the turbulence of the diffuse ISM.
We conclude that the decision of whether or not a density threshold is used has a significant and direct influence on the resulting VSFs. Indeed, it is a straightforward approach to focus the analysis on the actual area of interest. In observational studies such a threshold will anyway always be present as minimal collision rates for excitation or the sensitivity of detectors automatically introduce implicit density or intensity thresholds. Although we have only tested two specific setups in this context we have seen the significance of a proper choice of the density threshold, as well as a proper discussion of the obtained results considering the threshold as one of the defining parameters. We We note that future work could fruitfully also consider an upper density threshold to mimic the major effects of optical thickness even if true radiative transfer were probibitively expensive.
4.4 Density weighting
In this subsection, we discuss the effect of computing the VSF with or without including density weighting, relying on the results presented in Sect. 3.5. As long as the turbulence is dominated by the large scales, and a density threshold is used, considering the density weighting does not have a significant effect. However, as the clouds evolve the differences increase: the nonweighted VSFs never drop below 0.5. This is because the nonweighted VSF treats all cells equally, no matter whether the particular cell represents a dense element of the cloud centre or a diffuse element of the cloud’s edge, while the weighted VSF gives more weight to the matter within the smallscale, dense, collapsing regions. The kinetic energy is dominated by these regions. Thus, neglecting density weighting decouples the VSF from the kinetic energy distribution. A more exact treatment of this question is given by Kritsuk et al. (2013b) and Banerjee & Kritsuk (2017, 2018). This is particularly important at late times when smallscale collapse dominates.
Nevertheless, Fig. 3d illustrates that these differences do not prevent extended selfsimilarity from holding. Regardless of whether densityweighting is included, the values of Z remain similar, with similar responses to external driving, except for the features created when ζ(3) passes through zero in the densityweighted VSFs. This observation is true for all Jeans refinement levels, as Fig. 5 demonstrates.
Figure 5 summarises the comparison of ζ and Z measured with the densityweighted and nonweighted VSFs for all Jeans refinement levels (meaning the granularity used for modelling the turbulent motions of the gas, see Sect. 3.6 for more details). The figure clearly shows that the measurements only agree well for the highest refinement level with λ = 32Δx. However, we would need more data points to be sure that this correlation is indeed real. At lower refinement levels the measurements, as those used for the standard analysis and all other test scenarios but the one presented in Sect. 3.6, correlate less well with each other. The differences in the samples appear dominantly when the densityweighted ζ cease below ≈0.5, which is the global minimum for the nonweighted ζ. This means that none of the ζ computed in all clouds and refinement levels with the nonweighted VSF is measured to be below 0.5.
We conclude that deriving the VSF from smooth density distributions without considering densityweighting does not affect the behaviour ofζ and Z, as long as the turbulence is dominated by large scale flows, but it has a significant effect on the measurements when the small scales become dominant. The latter is particularly important as this finding has a directly impact on the conclusions drawn based on the scales and mechanisms that drive the turbulence based on the measured ζ. Not only does ζ become insensitive to the influence of gravitational contraction with time, the nonweighted VSFs also does not reflect when the majority of kinetic energy has been transferred to small scales. Furthermore, this emphasises the importance of taking optical depth effects into account, as single tracers covering limited density ranges may effectively provide statistics closer to the nonweighted VSFs.
4.5 Jeans length refinement
In Fig. 4 we see that the choice of refinement level has no significant influence on the measurements and evolution of both ζ and Z. The λ_{J} = 4Δx and λ_{J} = 8Δx models are in good agreement with each other. This means that, although refining Jeans lengths with 4 cells misses about 13% of kinetic energy, the effect on the structure and behaviour of the turbulence is rather small and not traced by the VSF analysis.
However, Fig. 4 shows that the agreement is rather poorer with λ_{J} = 32Δx, as the latter differs more from λ_{J} = 4Δx the higher the order of the VSF is. Following the explanations in Sect. 3.6, the behaviour of ζ and Z in the λ_{J} = 32Δx runs corresponds to the reaction of the cloud’s gas to a shock wave running through the cloud; caused by a SN that exploded before t = 0 Myr. Indeed one sees a SN at a distance of 172 pc at t = −1.11 Myr. As the power of the shock decreases rapidly with distance, the SN is too weak to effectively compress the gas within M3. This is why it was not detected in the less refined samples.
The SN explodes far below the midplane of the simulated disc galaxy, in a region without dense gas, so the blast wave remains strong as it propagates through the ISM. By the time the blast arrives at cloud M3, it is still energetic enough to impact the cloud with winds at velocities above 300 km s^{−1}, at the closer edge of the cloud. This causes an increase of VSFs at longer lag scales and the increase of ζ, as well as the drop in Z. However, in the less refined runs the arrival of the blast wave only contributes to the normal external driving for the gas’ turbulence. Only the λ_{J} = 32Δx runs can capture the fine structure of the shock driven into the cloud by the blast wave, which is why we detect this feature in these runs only. This emphasises the importance for studies like ours of properly resolving the full internal structure of realistic MCs.
We conclude that improving the resolution resolves details that can affect the VSF. However, the overall behaviour is already determined by our moderate resolution simulations.
Fig. 5
Comparison of ζ (top) and Z (bottom) measured based on densityweighted VSFs (abscissas) and nonweighted VSFs (ordinates). We note that the given values are based on M3 only. 
4.6 Comparison to observations
The majority of studies of VSFs in MCs are based on simulated data, as is the work presented in this paper. However, there are also some surveys that derive VSFs from observations, or whose data can be used to reconstruct the scaling properties of VSFs. In this section we discuss our results in the context of the observational studies that we list in Table 1.
Most of the velocity information derives from ^{12}CO and ^{13}CO observations of young starforming regions (e.g. Perseus and Taurus Padoan et al. 2003). We also consider observations of more evolved regions, such as those of the H II region NGC 6334 (Zernickel 2015), since the filaments in our simulated MCs fragment within the first 2 Myr Paper III, suggesting that our cloud may start to resemble such more evolved regions.
Figure 6 summarises the measured scaling exponents found in the literature, along with our fiducial results (Fig. 2a). We see that the observed values of ζ span within values of 0.24 to 1.18 across the first three orders. Furthermore, we see that the observed values are always positive, suggesting that none of the observed clouds show signs of being dominated by collapse, though that could be because the fastest flows lie in optically thick regions inaccessible to the observations.
Compared to these measurements, the distributions of our results show a large scatter across the parameter space. However, we also see that there is a significant fraction of values in our models that are consistent with the observational findings. These measurements belong roughly to the evolutionary stages of the modelled clouds after having evolved for roughly 1.0–3.9 Myr after the onset of selfgravity. Paper III finds that the clouds consist of a highly hierarchical structure that is dominated by already fragmenting filaments at this point. This means that the flows within the clouds experience a transition from cloudscale dominated, through filamentdominated, to corecollapse driven motions; which is exactly what we observe in the VSFs, as well. Consequently this suggests that the observed MCs, which show clear signs of embedded star formation activity, are in a similar stage where flows are dominated by the formation of hierarchical structures, (pre and proto)stellar cores or, in the case of NGC 6334, internal feedback.
We find that the interpretation of observational measurements is still difficult for several reasons:
 1.
We have already discussed in Sect. 4.2 that the transformation from 3D to 1D VSFs is not trivial, in particular when the studied flows are not isotropic. This is, for example, the case when the first structures (such as filaments or sheets) form, or the first cores collapse and accrete.
 2.
We have seen that interactions with SN shocks may trigger a preferred direction, that has the potential to strongly influence the measured 1D VSF. Although the influence of shock fronts on VSFs is transient compared to the lifetime of the entire MC, it is still long enough to mimic a quasisteady state in real observations. Observing typical shock tracers, such as SiO, may help to identify these situations.
 3.
We have neglected typical lineofsight effects that may have a significant influence on the measurements of the local standard of rest velocitywhose precision is crucial for this kind of study. Our projections ignore optical depth effects, and reflect velocities all the way through the clouds, including high column density regions of dynamical collapse where motions are fast at small scales. However, both ^{12}CO and ^{13}CO reach optical depth of unity at relatively low column densities. This means that the observed VSFs will only reflect the motions of the surface layers of dense MCs. Figure 2, as well as the figures in Appendix A demonstrate that neglecting certain regimes, in particular at small scales where column densities are highest, does have significant influence on the shape and scaling of the corresponding VSF. Therefore, singletracer observations are not suitable for studyingthe dynamical structure of MCs. For a proper VSF analysis it would be advisable to use a variety of tracers tocover the different phases of the clouds, as well as to populate the statistics of lag distances more completely.
 4.
Only a small fraction of the listed observational studies in Table 1 aimed to measure the VSFs of the respective objects directly. In the majority of cases, the focus of the investigations was on the general budget of kinetic energy within the MC, as well as the question whether those clouds follow Larson’s sizevelocity relation (Eq. (1)). It is unclear whether the difference between a relation of the lag distance of two particles and their relative lineofsight velocity and the connection between the size of the entire MC and the velocity dispersion of the contained gas has always been considered.
We recommend that both theorists and observers discuss in more detail how observational studies may use VSFs in the future. From the theoretical pointofview, full line radiative transfer calculations are required to better evaluate observational biases and simple projection effects. This requires observations with a high spatial resolution of the respective MC for a wide range of lag scales and good statistics for fitting the scaling of VSF, as well as lines with welldefined lineofsight velocities, ideally, optically thin lines of intermediate and highdensity tracers.
Summary of observed ζ and Z in the literature.
Fig. 6
Time evolution of scaling exponent ζ(p) of the pth order VSF. The panels show the measured scaling exponents for the first to third order (left to right) for M3 (red), M4 (green), M8 (blue). The shaded area within all panels marks the range of observed values at the respective order of VSF (see Table 1). 
5 Summary and conclusions
In this paper, we analyse the VSFs of MCs that have formed within 3D magnetohydrodynamical, adaptive mesh refinement, FLASH simulations of the selfgravitating, SNdriven ISM by Paper II, including both density weighting and a density cutoff. The main results are as follows.
 1.
The scaling of the densityweighted VSFs depends on a combination of turbulence and more coherent processes such as SN blast wave impacts and gravitational contraction. We find that the powerlaw scaling ζ of 3D densityweighted VSFs reflects the development of gravitational contraction, while the extended selfsimilarity scaling Z reveals interactions of clouds with largescale blast waves.
 2.
The two different proposed explanations for Larson’s sizevelocity relationship, a turbulent cascade and gravitational contraction, appear to apply to different stages in the evolution of MCs, as well as different observational techniques. It appears coincidental that they have the same functional relationship of length to velocity, which has led to confusion of one with the other.
MCs dominated by uniform turbulence show a firstorder VSF with ζ(1) ≃ 1∕2. The same result can be found for clouds undergoing strong gravitational contraction by computing the VSF without density weighting (or, most likely, in the presence of optical depth effects), which is dominated by the lowdensity, turbulencedominated outer regions of clouds. At the initial time, though, we measure values of ζ(2) = 1 or larger in the densityweighted VSFs, that are inconsistent with turbulence models or simulations that predict ζ(2) = 0.74 (Eq. (6)). However, these initial values are affected by numerical dissipation and tend to decrease as the resolution in dense regions increases in the zoominruns.
Examining the overall velocity dispersion of gravitationally dominated clouds undergoing star formation, on the other hand, reflects the dynamics of gravitational collapse. In this case, the cloud shows a shallow or even inverted VSF dependence ζ(1) ≲ 0. This reflects strong flows at small scales. However, such gravitationally contracting clouds were shown by Paper I to have an overall squareroot velocityradius relationship (Eq. (1)) given by freefall or virial equilibrium (which differ by only 2^{1∕2}, as notedby BallesterosParedes et al. 2011b).
 3.
As long as the MC is not affected by a shock, Z agrees well with predicted values for supersonic flows, even as gravitational collapse proceeds.
 4.
We test the influence of Jeans refinement on the VSFs. We find that the absolute amount of kinetic energy does not influencethe evolution of ζ and Z, but thatbetter resolution of external shocks can produce changes in both quantities.
 5.
Comparison of 3D to 1D VSFs shows differences in detail, but qualitative agreement in the behaviour of both ζ and Z, in particular when gravity dominates gas dynamics. Thus, observed 1D VSFs can be useful diagnostics in gravitationally bound and contracting regions. On the other hand, differences arise when strong transverse flows or shocks dominate the velocity field.
 6.
We calculate cloud VSFs using a density threshold to isolate the cloud material, as would characteristically happen in an observation of molecular material. Without such a threshold, our VSFs are dominated by the diffuse ISM. In that case, the extended selfsimilarity scaling Z lies just below the value predicted for isotropic, incompressible turbulence by She & Lévêque (1994). This is consistent with the low Mach number in the hot, diffuse, ISM filling most of the volume of our simulation. Yet, the actual values of ζ in the low density gas are about four times higher than those predicted by both She & Lévêque (1994) and Boldyrev (2002) due to some combination of numerical dissipation and the multiphase nature of the medium, which reduces velocities at small scale in dense regions (see Sect. 4.3).
 7.
We investigate the influence of defining the VSF with and without density weighting. We find that the qualitative behaviour is traced by both approaches. However, the scaling of the nonweighted VSF ζ is always positive, not falling nearly as far as for the densityweighted VSF. The densityweighted VSF reflects the kinetic energy distribution better as gravitational collapse proceeds to smaller and smaller scales. (We note that in, for example, CO observations, optical depth effects may obscure this behaviour.)
 8.
We compare our results with measurements of both ζ and Z in observational studies. We see that our findings are generally consistent with with observations within periods during which the clouds’ flows are influenced by both turbulent flows and global gravitational contraction, including strong structure formation and starting fragmentation. This reflects the conditions of embedded star formation activity within observed MCs.
Our analysis shows that VSFs are useful tools for examining the driving source of turbulence within MCs. However, studies that use VSFs need to precisely review the assumptions and parameters included in their analysis as these can have a significant influence on the results.
For our simulated clouds, the VSFs illustrate that gravitational contraction dominates the evolution of the clouds. During contraction, the VSF scaling exponent ζ(p) drops in value and can even become negative as kinetic energy concentrates on small scales. Nevertheless, the extended selfsimilarity scaling parameters Z(p) continue to agree with the analytic prediction for compressible turbulence except for short periods during which SN blast waves increase power on multiple scales. Because such blast waves are neither homogeneous nor isotropic, they often lead to transient nonpower law scaling of the VSFs, and thus strong departures from uniform turbulent behaviour of Z(p).
Acknowledgements
We thank the anonymous referee for a detailed report that led to an improved and better focussed exposition, particularly in the discussion of Larson’s relation, as well as the inclusion of the Appendices. M.M.M.L. received support from US NSF grants AST1109395 and AST1815461, and thanks the A. von HumboldtStiftung for support. J.C.I.M. was additionally supported by the Deutsche Forschungsgemeinschaft (DFG) via the Collaborative Research Center SFB 956 “Conditions and Impact of Star Formation” (subproject C5) and the DFG Priority Program 1573 “The physics of the interstellar medium”.
Appendix A Computation and fitting procedures
In this section we provide more details on how we compute the VSFs and their scaling parameters. As described in Sect. 2.2 the discussions in the main part of this manuscript are based on VSFs that are computed from average relative velocities. This means the following:
We mapthe 3D FLASH adaptive mesh refinement data of the original simulations (Papers I; II) onto uniform grid data cubes. The cubes simulate volumes of 40 pc on a side with 0.1 pc zones centred on the centres of the molecular clouds.
The FLASH simulations take ten levels of resolution into account; and the decision which region is resolved up to which level depends on the density of the respective grid cell (Paper II). Therefore, the original simulations cover a range of resolutions from 30.4 pc in regions of z > 10 kpc above and below the galactic midplane down to 0.06–0.10 pc within (100 pc)^{3} boxes around the centre of mass of the three molecular clouds presented in this manuscript, as well as in Paper II and Paper III. Therefore, the resolution we have chosen for the data cubes corresponds to either the highest resolved level offered by the simulations (in the case of M8) or a slight coarsening of the most resolved level (in the cases of M3 and M4). The zoomin regions modelled at the highest resolution level extend well beyond the volume we use in our analysis. We have adopted this size as it is sufficient for covering theclouds’ volumes, and thus more efficient for further analyses. In addition, this choice protects us from resolution edge effects as the transitions from this resolution level to the next lowest are several tens of zones away.
No matter whether a density threshold is applied or not, the (40 pc)^{3} cubes still include too much data to compute complete velocity structure functions including all lags from all points. Therefore, to derive the VSFs we coarsen the grid of projected lag distances, ℓ_{i} = ℓ_{i} to cover the range 0.8–30 pc with only 40 equidistant bins; relative to each of the starting coordinates of our samples.
After mapping the data onto the uniform grid, we apply two approaches: The first considers only zones above the density threshold, n_{cloud}. These zones represent the starting points x (see Eqs. (2)–(8)). Our routine calculates the lag distances between these zones and every other zone in the sample, as well as the relative velocities of the gas within the given zones. The individual lag distances are binned using spherical shells around the starting zones that range from inner radii ℓ_{i} to outer radii ℓ_{i+1}. By doing so, we compute the discrete VSFs presented in the main part of this paper with the relative velocities and product of densities, ρ(x) ⋅ ρ(x + ℓ), measured from zones within the individual shells.
Fig. A.1
Data presented in Fig. 2a, with shaded areas behind the data representing the error ranges of the computed ζ (top) and Z (bottom). 
The second approach targets the case when we do not apply any density threshold (i.e. setting n_{cloud} = 0). In this case, we use a random number generator (random.rand from numpy) to choose 5% of the total zones, and do the analysis on them. We emphasise that this does not mean that we only calculate the relative velocities between these zones. Rather this subsample of zones represent the starting vectors x to which the velocities of all other zones x + ℓ in the same cube are compared to. This way we reduce the risk of ignoring or emphasising any spatial direction or angle. As it is too computationally expensive to derive all relative velocities between all zones within discrete shells, as we have done in the first approach, we derived a discrete distribution of relative velocities as function of lag distance using a discrete fast Fourier transform (FFT). These distributions are based on the same grid of lag distances we have already utilised for the spherical shells above. Therefore, we can use the results of the FFT in the same way as the results from the first approach to derive the VSFs.
With both approaches we obtain a set of discrete descriptions of VSFs as a function of lag and time (by computing VSFs for successive code outputs). In order to derive the scaling parameters ζ of the VSFs as a function of time and order, we fit the powerlaw relation presented in Eq. (10) to the measured VSFs, using the python curve_fit package. Due to the rather irregular behaviour of our VSFs at larger scales, we define the weighting function as (A.1)
Although the average radii of our clouds are, on average, larger than 8 pc we choose this limit due to the variable behaviour of the VSFs at scales of that size and larger. Investigation of different weightings would be a fruitful topic of further investigation.
As can be seen in Fig. 1, as well as in the figures shown in Appendix B, the shape of VSFs changes over time as different forces act, as we explain in Sect. 3. The similarity parameters, Z, are computed by applying Eq. (9) on the results of the fitting procedure, and results are presented in Sect. 3, as well as in Appendix B.
The fit also provides the χ^{2} errors for the measured values of ζ. In Fig. A.1 we show a reduced version of Fig. 2(a), where we only plot the time evolution of ζ for all three clouds, along with shades of the same colours of the respective lines that represent the errors. We see that the relative errors, Δζ∕ζ mostly remain within a range of 5–12%. The errors of Z are computed by Gaussian error propagation
In general, the relative errors of Z are, as well, around 10%, though we do see exceptions with very large errors. The reason for these is that at these times ζ(3) approaches zero, causing Eq. (A.3) to diverge.
Appendix B Inertial range
For our analysis it is essential to verify that there is a reasonable range of scales ℓ below the driving range within which the simulations are not dominated by numerical effects, so that they resolve the inertial range of any turbulent cascade. We note that the fitting region for our VSFs starts at lags greater than eight zones, ensuring that the numerical dissipation range lies at smaller scales. To understand the nonpower law behaviour that we nonetheless find, in the following subsections we offer more examples that show VSFs of all three clouds at different times and considering our different analysis approaches. We focus on the standard analysis in Appendix B.1, on the analysis neglecting the density threshold in Appendix B.2, and on the impact of varying the resolution in collapsing regions in Appendix B.3.
B.1 Standard analysis
Figure B.1 extends the data presented in Sect. 3.2 and discussed in Sect. 4.1. The figure uses the same format as Fig. 1. The straight lines within the plots indicate the powerlaw relation that we have fitted onto the VSFs, considering the range 0.8 ≤ ℓ ≤ 8 pc.
We see that in most of the cases the VSFs are in good agreement with the described powerlaw relation within the fitted ranges (e.g. Fig. 1). However, there are cases when the VSF is not well reproduced by a simple, single powerlaw function, such as M3 at t = 3.0 Myr in Fig. 1. This appears to occur particularly when the dominant mechanism driving turbulence changes from largescale driving to internal contraction. At this point the larger scales of the clouds (larger ℓ) are still dominated by external driving, while the smaller scales (smaller ℓ) start to accelerate mostly due to gravitational fragmentation and infall motions. The consequence is that the actual VSF is a superposition of two processes that amplify the relative motions of the gas differently and on different scales. A single powerlaw cannot describe this scenario adequately.
B.2 No density threshold
Figures B.2 and B.3 extend the set of VSFs for the box without a density threshold presented in Sect. 3.4 and discussed in Sect. 4.3. Here we also find problems describing the behaviour of smallscale motions with ℓ ≲ 2 pc. In this case we primarily capture the turbulent motions within the lowdensity ISM. Contrary to the modelled MCs, the ISM is not organised in hierarchical structures, and the turbulence within the ISM is predominantly driven by SN explosions. These produce structured correlations at the large scales of entire blast waves as well as the small scales across shockfronts. Yet, we see that we can fit the VSFs with a power law within the intermediate lag scale ranges (2 pc ≲ ℓ ≲ 8 pc), suggesting this range of length scales is dominated by a turbulent cascade.
B.3 Varying refinement
Figures B.4 and B.5 extend the set of VSFs presented for clouds with varying amounts of refinement in gravitationally collapsing regions in Sect. 3.6 and discussed in Sect. 4.5. A final example of VSFs not following a power law at intermediate scales is given by cloud M3 at t = 0.8 Myr and with λ_{J} = 32Δx in Fig. B.4. This represents the VSF of a cloud that is interacting with a SN blast wave. In this case the maximal amplification is neither at the scale of the cloud, where the turbulence is driven by external sources, nor on small scales where gravitational contraction acts. Instead the local maximum is the intermediate scale. Considering the morphology of the cloud and the cloud’s environment this can only occur when the shock front of a SN is currently propagating through the cloud. Thus, the VSF here is a superposition of three driving mechanisms: first, the external largescale driving; second, selfgravity that leads to contractive motions; and finally the shock jump that injects kinetic energy as it moves through the cloud. The effect of the shock, however, is only local and shortlived as the injected turbulence decays quickly (compare with M3 at t = 1.2 Myr and with λ_{J} = 32Δx in Fig. B.4).
Fig. B.4
Additional examples of VSFs, based on data of M3 with density threshold, at different refinement levels (left to right) λ = 4Δx, λ = 8Δx, and λ = 8Δx as function of lag scale ℓ and order p. The examples are given for three different time steps, namely (top to bottom) t = 0.0, 0.8, and 1.2 Myr. The dots (connected by dashed lines) show the values computed from the simulations. The solid lines represent the powerlaw relation fitted to the respective structure functions. 
References
 BallesterosParedes, J., Hartmann, L. W., VázquezSemadeni, E., Heitsch, F., & ZamoraAvilés, M. A. 2011a, MNRAS, 411, 65 [NASA ADS] [CrossRef] [Google Scholar]
 BallesterosParedes, J., VázquezSemadeni, E., Gazol, A., et al. 2011b, MNRAS, 416, 1436 [Google Scholar]
 Banerjee, S., & Galtier, S. 2013, Phys. Rev. E, 87, 013019 [NASA ADS] [CrossRef] [Google Scholar]
 Banerjee, S., & Kritsuk, A. G. 2017, Phys. Rev. E, 96, 053116 [NASA ADS] [CrossRef] [Google Scholar]
 Banerjee, S., & Kritsuk, A. G. 2018, Phys. Rev. E, 97, 023107 [NASA ADS] [CrossRef] [Google Scholar]
 Benzi, R., Ciliberto, S., Tripiccione, R., et al. 1993, Phys. Rev. E, 48, R29 [NASA ADS] [CrossRef] [Google Scholar]
 Benzi, R., Biferale, L., Fisher, R., Lamb, D. Q., & Toschi, F. 2010, J. Fluid Mech., 653, 221 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Boldyrev, S. 2002, ApJ, 569, 841 [NASA ADS] [CrossRef] [Google Scholar]
 Boneberg, D. M., Dale, J. E., Girichidis, P., & Ercolano, B. 2015, MNRAS, 447, 1341 [NASA ADS] [CrossRef] [Google Scholar]
 Brunt, C. M., & Heyer, M. H. 2013, MNRAS, 433, 117 [NASA ADS] [CrossRef] [Google Scholar]
 Brunt, C. M., Heyer, M. H., & Mac Low, M.M. 2009, A&A, 504, 883 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Burkhart, B., Collins, D. C., & Lazarian, A. 2015, ApJ, 808, 48 [NASA ADS] [CrossRef] [Google Scholar]
 Chira, R.A., IbáñezMejía, J. C., Mac Low, M.M., & Henning, T. 2018a, How do Velocity Structure Functions Trace Gas Dynamics in Simulated Molecular Clouds? (New York: American Museum of Natural History Research Library) [Google Scholar]
 Chira, R.A., Kainulainen, J., IbáñezMejía, J. C., Henning, T., & Mac Low, M.M. 2018b, A&A, 610, A62 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dekel, A., & Krumholz, M. R. 2013, MNRAS, 432, 455 [NASA ADS] [CrossRef] [Google Scholar]
 Elmegreen, B. G. 1993, in Protostars and Planets III, eds. E. H. Levy & J. I. Lunine (Tucson: University of Arizona Press), 97 [Google Scholar]
 Elmegreen, B. G. 2007, ApJ, 668, 1064 [NASA ADS] [CrossRef] [Google Scholar]
 Falgarone, E., Pety, J., & HilyBlant, P. 2009, A&A, 507, 355 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Federrath, C., Sur, S., Schleicher, D. R. G., Banerjee, R., & Klessen, R. S. 2011, ApJ, 731, 62 [NASA ADS] [CrossRef] [Google Scholar]
 Fleck, Jr. R. C. 1980, ApJ, 242, 1019 [NASA ADS] [CrossRef] [Google Scholar]
 Fryxell, B., Olson, K., Ricker, P., et al. 2000, ApJS, 131, 273 [NASA ADS] [CrossRef] [Google Scholar]
 Galtier, S., & Banerjee, S. 2011, Phys. Rev. Lett., 107, 134501 [NASA ADS] [CrossRef] [Google Scholar]
 Gnedin, O. 2015, IAU General Assembly, 22, 2256326 [NASA ADS] [Google Scholar]
 Goodman, A. A., Barranco, J. A., Wilner, D. J., & Heyer, M. H. 1998, ApJ, 504, 223 [NASA ADS] [CrossRef] [Google Scholar]
 Gotoh, T., Fukayama, D., & Nakano, T. 2002, Phys. Fluids, 14, 1065 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Hartmann, L., BallesterosParedes, J., & Heitsch, F. 2012, MNRAS, 420, 1457 [Google Scholar]
 Heyer, M. H., & Brunt, C. M. 2004, ApJ, 615, L45 [NASA ADS] [CrossRef] [Google Scholar]
 Heyer, M. H., & Brunt, C. 2007, IAU Symp., 237, 9 [NASA ADS] [Google Scholar]
 Heyer, M., & Dame, T. M. 2015, ARA&A, 53, 583 [NASA ADS] [CrossRef] [Google Scholar]
 Heyer, M., Krawczyk, C., Duval, J., & Jackson, J. M. 2009, ApJ, 699, 1092 [NASA ADS] [CrossRef] [Google Scholar]
 IbáñezMejía, J. C., Mac Low, M.M., Klessen, R. S., & Baczynski, C. 2016, ApJ, 824, 41 [NASA ADS] [CrossRef] [Google Scholar]
 IbáñezMejía, J. C., Mac Low, M.M., Klessen, R. S., & Baczynski, C. 2017, ApJ, 850, 62 [NASA ADS] [CrossRef] [Google Scholar]
 Kolmogorov, A. 1941, Akademiia Nauk SSSR Doklady, 30, 301 [Google Scholar]
 Kritsuk, A. G., Lee, C. T., & Norman, M. L. 2013a, MNRAS, 436, 3247 [NASA ADS] [CrossRef] [Google Scholar]
 Kritsuk, A. G., Wagner, R., & Norman, M. L. 2013b, J. Fluid Mech., 729, R1 [CrossRef] [Google Scholar]
 Kritsuk, A. G., Wagner, R., & Norman, M. L. 2015, ASP Conf. Ser., 498, 16 [NASA ADS] [Google Scholar]
 Krumholz, M. R., Bate, M. R., Arce, H. G., et al. 2014, Protostars and Planets VI (Tucson: University of Arizona Press), 243 [Google Scholar]
 Larson, R. B. 1981, MNRAS, 194, 809 [NASA ADS] [CrossRef] [Google Scholar]
 Mac Low, M.M. 2003, in Turbulence and Magnetic Fields in Astrophysics, Lect. Notes Phys., (Berlin: Springer Verlag), eds. E. Falgarone & T. Passot, 614, 182 [NASA ADS] [CrossRef] [Google Scholar]
 Mac Low, M.M., & Klessen, R. S. 2004, Rev. Mod. Phys., 76, 125 [NASA ADS] [CrossRef] [Google Scholar]
 McKee, C. F., & Zweibel, E. G. 1992, ApJ, 399, 551 [NASA ADS] [CrossRef] [Google Scholar]
 Miesch, M. S., & Bally, J. 1994, ApJ, 429, 645 [NASA ADS] [CrossRef] [Google Scholar]
 Miyamoto, Y., Nakai, N., & Kuno, N. 2014, PASJ, 66, 36 [NASA ADS] [CrossRef] [Google Scholar]
 Padoan, P., Boldyrev, S., Langer, W., & Nordlund, Å. 2003, ApJ, 583, 308 [NASA ADS] [CrossRef] [Google Scholar]
 Padoan, P., Juvela, M., Kritsuk, A., & Norman, M. L. 2006, ApJ, 653, L125 [NASA ADS] [CrossRef] [Google Scholar]
 Padoan, P., Pan, L., Haugbølle, T., & Nordlund, Å. 2016, ApJ, 822, 11 [NASA ADS] [CrossRef] [Google Scholar]
 RomanDuval, J., Federrath, C., Brunt, C., et al. 2011, ApJ, 740, 120 [NASA ADS] [CrossRef] [Google Scholar]
 Schmidt, W., Federrath, C., & Klessen, R. 2008, Phys. Rev. Lett., 101, 194505 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Seifried, D., Walch, S., Girichidis, P., et al. 2017, MNRAS, 472, 4794 [Google Scholar]
 She, Z.S., & Lévêque, E. 1994, Phys. Rev. Lett., 72, 336 [Google Scholar]
 Solomon, P. M., Rivolo, A. R., Barrett, J., & Yahil, A. 1987, ApJ, 319, 730 [NASA ADS] [CrossRef] [Google Scholar]
 Tan, J. C., Shaske, S. N., & Van Loo, S. 2013, IAU Symp., 292, 19 [NASA ADS] [Google Scholar]
 Truelove, J. K., Klein, R. I., McKee, C. F., et al. 1998, ApJ, 495, 821 [NASA ADS] [CrossRef] [Google Scholar]
 Turk, M. J., Oishi, J. S., Abel, T., & Bryan, G. L. 2012, ApJ, 745, 154 [NASA ADS] [CrossRef] [Google Scholar]
 VázquezSemadeni, E., Ryu, D., Passot, T., González, R. F., & Gazol, A. 2006, ApJ, 643, 245 [NASA ADS] [CrossRef] [Google Scholar]
 Zernickel, A. 2015, PhD Thesis, Köln, Germany [Google Scholar]
All Tables
All Figures
Fig. 1
Examples of VSFs from models (left to right) M3, M4, and M8 as function of lag scale ℓ and order p, based on data with density threshold. The examples are given for times (top to bottom) t = 1.0 Myr, 3.0 Myr, and 4.2 Myr. The dots (connected by dashed lines) show the values computed from the simulations. The solid lines represent the powerlaw relation fitted to the respective structure functions. 

In the text 
Fig. 2
Time evolution of scaling exponent ζ(p) of the pth order VSF. Panels a–d: the measurements for M3 (left), M4 (middle), and M8 (right), respectively. Of these, panel a represents the standard analysis while the other rows illustrate the results of different variations in the analysis as noted in the figures. Panel e: the values of ζ measured within M3 as a function of the Jeans refinement level the cloud has been modelled with. We note that these more expensive runs were not run for as long as the fiducial run. In all panels, the grey dotted vertical lines mark the times than a SN explodes within 100 pc of the cloud, while the blue areas indicate periods of enhanced mass accretion onto the clouds. The coloured horizontal lines show the predicted values for incompressible turbulence (dashdotted lines; She & Lévêque 1994) and for highly compressible, supersonic turbulence (dashed lines; Boldyrev 2002). 

In the text 
Fig. 3
Like Fig. 2, but for the measured selfsimilarity parameter Z = ζ(p)∕ζ(3) of the pth order VSF. 

In the text 
Fig. 4
Comparison of VSF scaling exponents, ζ (left), and selfsimilarity parameters, Z (right), depending on the Jeans refinement of the simulation runs the data are based on. The abscissas give values from λ_{J} = 4Δx, while the ordinates give values from λ_{J} = 8Δx (top) and λ_{J} = 32Δx (bottom). All data points refer to the M3 cloud and represent different lags in the same time step in the respective simulations. 

In the text 
Fig. 5
Comparison of ζ (top) and Z (bottom) measured based on densityweighted VSFs (abscissas) and nonweighted VSFs (ordinates). We note that the given values are based on M3 only. 

In the text 
Fig. 6
Time evolution of scaling exponent ζ(p) of the pth order VSF. The panels show the measured scaling exponents for the first to third order (left to right) for M3 (red), M4 (green), M8 (blue). The shaded area within all panels marks the range of observed values at the respective order of VSF (see Table 1). 

In the text 
Fig. A.1
Data presented in Fig. 2a, with shaded areas behind the data representing the error ranges of the computed ζ (top) and Z (bottom). 

In the text 
Fig. B.1
As Fig. 1, but plotting relation between S_{ℓ}/ℓ as function of lag scale l and order p. 

In the text 
Fig. B.2
As Fig. 1, but based on data without density threshold. 

In the text 
Fig. B.3
As Fig. B.1, but based on data without density threshold. 

In the text 
Fig. B.4
Additional examples of VSFs, based on data of M3 with density threshold, at different refinement levels (left to right) λ = 4Δx, λ = 8Δx, and λ = 8Δx as function of lag scale ℓ and order p. The examples are given for three different time steps, namely (top to bottom) t = 0.0, 0.8, and 1.2 Myr. The dots (connected by dashed lines) show the values computed from the simulations. The solid lines represent the powerlaw relation fitted to the respective structure functions. 

In the text 
Fig. B.5
As Fig. B.4, but plotting relation between S_{ℓ}/ℓ as function of lag scale l and order p. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.