Supernova Driving. I. The Origin of Molecular Cloud Turbulence

Turbulence is ubiquitous in molecular clouds (MCs), but its origin is still unclear because MCs are usually assumed to live longer than the turbulence dissipation time. Interstellar medium (ISM) turbulence is likely driven by SN explosions, but it has never been demonstrated that SN explosions can establish and maintain a turbulent cascade inside MCs consistent with the observations. In this work, we carry out a simulation of SN-driven turbulence in a volume of (250 pc)$^3$, specifically designed to test if SN driving alone can be responsible for the observed turbulence inside MCs. We find that SN driving establishes a velocity scaling consistent with the usual scaling laws of supersonic turbulence, suggesting that previous idealized simulations of MC turbulence, driven with a random, large-scale volume force, were correctly adopted as appropriate models for MC turbulence, despite the artificial driving. We also find that the same scaling laws extend to the interior of MCs, and that the velocity-size relation of the MCs selected from our simulation is consistent with that of MCs from the Outer-Galaxy Survey, the largest MC sample available. The mass-size relation and the mass and size probability distributions also compare successfully with those of the Outer Galaxy Survey. Finally, we show that MC turbulence is super-Alfv\'{e}nic with respect to both the mean and rms magnetic-field strength. We conclude that MC structure and dynamics are the natural result of SN-driven turbulence.


INTRODUCTION
Understanding molecular cloud (MC) turbulence is key to understanding the star formation process, because supersonic turbulence is ubiquitous in MCs and drives their fragmentation into stars. Supersonic turbulence has been studied extensively in the context of star formation, and its statistical properties are at the core of recent models of the star formation rate (Krumholz & McKee 2005;Hennebelle & Chabrier 2011; and the stellar initial mass function Hennebelle & Chabrier 2008;Hopkins 2012). Most numerical studies of supersonic turbulence have used a random large-scale force to drive the turbulence, as customary in the turbulence literature. It remains to be shown that such an idealized external force is a good approximation of the actual large-scale processes driving the interstellar-medium (ISM) turbulence. While the simulations use a volume force that penetrates the interior of the fluid, as long-range forces do (e.g. gravity and magnetic fields), the real ISM driving forces may be surface forces, such as large-scale shocks from spiral arms or su-pernova (SN) bubbles. The effect of different types of surface forces, or of a combination of volume and surface forces on the turbulence and, thus, on star formation, has not been studied systematically.
The turbulence in MCs is also key to understanding their origin. The generation and maintenance of MC turbulence must be an integral part of the cloud formation process, because most of the energy of an MC is in the form of turbulent kinetic energy. For example, the great majority of small and intermediate-mass MCs are known to have rather large virial parameters (see §10 andHeyer et al. 2001, 2009), thus their turbulent energy is large enough to form and disperse them in a few dynamical times. The same may be true also for the most massive MCs, even if their virial parameter tends to be of order unity. The spatial and velocity structures of MCs follow power laws that span all scales from the smallest to the most massive clouds, suggesting a universal origin of clouds and cloud turbulence.
In this work, we adopt the viewpoint that MC turbulence is just one component of the general ISM turbulence. The question of the origin of MC turbulence is thus turned into the more general question of how the ISM turbulent energy is shared among its different gas phases. If we demonstrate that the total large-scale turbulent energy of the ISM is dynamically consistent with the turbulent energy in its dense phase, no extra energy source specific to MCs is needed.
It is generally accepted that SN explosions dominate the energy budget of star-forming galaxies at MC scales, although large-scale gravitational instabilities in galactic disks (e.g. Elmegreen et al. 2003;Bournaud et al. 2010) and gas compression in spiral density waves (e.g. Semenov et al. 2015) may also contribute to the turbulence. The Kennicut-Schmidt star formation law of disk galaxies, despite giving a large gas-consumption timescale of order 1 Gyr, corresponds to an energy input from SN explosions that exceeds the turbulence dissipation rate of the ISM of those galaxies. The analysis of HI maps of nearby face-on galaxies also leads to the conclusion that SN feedback is responsible for the observed HI line width, except in the disk outskirts (Tamburro et al. 2009;Stilp et al. 2013;Ianjamasimanana et al. 2015). Detailed modeling of the disk vertical balance has also shown that SN feedback can maintain the ISM turbulence that determines the disk scale height, resulting in self-regulated star formation (Ostriker et al. 2010;Ostriker & Shetty 2011;Faucher-Giguère et al. 2013), a picture that may apply also to galaxies with high star formation rate at redshift z = 1 − 3 (Lehnert et al. 2013).
Prior to these studies, galactic-fountain simulations had already demonstrated that SN explosions can drive the observed ISM turbulence (de Avillez & Breitschwerdt 2005;Joung et al. 2009). However, these simulations barely resolved the formation of the most massive clouds, and could not resolve their internal structure and dynamics in order to compare with observed MC properties. Furthermore, the evolution of individual SN explosions was not resolved with high-enough spatial resolution to study their interaction with individual MCs.
Recent numerical works have studied the momentum injection by individual SN explosions into the ambient medium, assuming realistic ISM density and temperature fluctuations Martizzi et al. 2015;Iffrig & Hennebelle 2015;Kim & Ostriker 2015). These studies are useful for the derivation of feedback models to be implemented as sub-grid physics in galaxy formation simulations, but do not address the problem of the generation and maintenance of MC turbulence by SN explosions.
In this work, we focus on a smaller-scale region than in galactic-fountain simulations, and use high-enough numerical resolution to model the MCs, their internal structure and dynamics, and the evolution of individual SN bubbles with sub-pc resolution (see §2). We also select clouds from the simulation, as connected regions above a threshold density ( §3), to study their turbulence and to compare with observed MCs. First we study statistical properties over the whole computational volume, such as total energies ( §4), power spectra ( §5) and velocity structure ( §6.1). Then we present properties of individual clouds selected from the simulation, such as the velocity structure ( §6.2), the virial parameter ( §7), the cloud lifetime ( §8) and the magnetic field strength ( §9). Finally, we compare the properties of the clouds selected from the simulation with those of MCs from the Five College Radio Astronomy Observatory (FCRAO) Outer Galaxy Survey ( §10).

NUMERICAL SETUP
The simulation is carried out with the Ramses AMR code (Teyssier 2002). We refer to Padoan et al. (2014) for a brief description of our version of Ramses. What is new here, relative to the simulation in Padoan et al. (2014), is the use of the full energy equation (instead of assuming an isothermal equation of state), and the inclusion of SN explosions and tracer particles, besides the much larger physical size of the computational volume.
We simulate a cubic region of size L box = 250 pc (large enough to contain a few turbulence correlation lengths, while small enough to allow sub-pc resolution), with a mean density of 5 cm −3 (corresponding to a column density of 30 M pc −2 and a total mass of 1.9 × 10 6 M ) and a Type II SN rate of 6.25 Myr −1 (or a galactic rate of 100 Myr −1 kpc −2 , if all SN explosions occurred within the vertical extent of our box). We do not consider Type Ia SNe, because of their lower rate and higher scale height, and because we distribute SN explosions randomly, so our SN rate could also be interpreted as the sum of the Type II and Type Ia rates.
These rate and column density values are consistent with the Kennicutt-Schmidt relation, and our computational volume may be viewed as a dense section of a spiral arm. For example, the total column density in the Perseus arm of the Milky Way is 23 M pc −2 (Heyer & Terebey 1998). However, we do not include a galactic gravitational field, and adopt periodic boundary conditions in all directions, so vertical stratification and outflows of hot gas are neglected. We have chosen this idealized setup because one of the motivations of this work is to relate the statistical properties of SN-driven turbulence to previous studies of randomly-driven, supersonic, isothermal turbulence that were carried out on periodic boxes without stratification. Furthermore, we relate the velocity scaling to theoretical predictions that also neglect complications such as gravity and stratification. 1 Besides the pdV work, and the thermal energy introduced to model SN explosions, our energy equation adopts uniform photoelectric heating as in Wolfire et al. (1995), with efficiency = 0.05 and the FUV radiation field of Habing (1968) with coefficient G 0 = 0.6, chosen to obtain temperature distributions consistent with those from the comprehensive simulations by . Because the code conserves total energy, kinetic and magnetic dissipations are included self-consistently as energy sources (this dissipation is purely numerical, as we do not include viscosity or resistivity explicitly). We use a tabulated optically thin cooling function constructed from the extensive compilation by Gnedin & Hollon (2012), based on 75 million runs of the Cloudy code (Ferland et al. 1998) to sample a large range of conditions, and from which the results have been made publicly available as a Fortran code with accompanying database. All relevant atomic transitions are included in the Cloudy runs. Although available in Cloudy, molecu-1 Star-formation simulations with gravity have shown that, under reasonable conditions (e.g. MCs not collapsing as a whole), gravity does not affect the velocity scaling (Federrath & Klessen 2013) Fig. 1.-Phase diagram of gas pressure versus density based on 11 snapshots covering uniformly the time period when gravity is included in the simulations, t = 45 to 56 Myr. The gray scale represents the square root of the volume fraction at a given pressure and density. lar cooling is not included because the runs are restricted to a single computational zone, with a negligible column density, to enforce the optically thin case. Above a temperature of 100 K, atomic cooling is dominant up to densities of 10 6 cm −3 (e.g. Neufeld et al. 1995). At lower temperature and high densities, molecular cooling should be included. However, molecular cooling and cosmic-ray heating are neglected, and their thermal balance in very dense gas is emulated by clamping the resulting drop in temperature at 10 K. As pointed out by Gnedin & Hollon (2012), including the balance of molecular cooling and cosmic-ray heating in such a treatment does not make much physical sense, since the balance of these processes at high densities crucially depends on radiative transfer effects; in particular the absorption of UV radiation by small-scale high-density cloud structures. In the absence of radiative transfer, and given the optically thin assumption, we approximate the UV shielding in MCs by tapering off the photoelectric heating exponentially above a number density of 200 cm −3 (assuming a characteristic size of 1 pc for MC structures at our critical density, corresponding to a critical visual extinction of 0.3 mag -Franco & Cox 1986). Figure 1 shows the phase diagram of gas pressure versus density sampled over the last 11 Myr of the simulation. The horizontal feature around densities of a few 10 −22 g cm −3 is a consequence of the approximation of self-shielding. In a real MC, or a model where the absorption of UV-radiation by the filamentary structure of dense gas and dust is taken into account more realistically, the transition to opaqueness would take place at different densities at different locations, and similar local phase diagram features would be washed out. Comparing our current phase diagram with the ones in , we actually find the largest discrepancy not to be near that horizontal feature, but rather at higher densities where, for some reason, their balance of heating and cooling results in a temperature of about 30 K, instead of our assumed value of 10 K.
The simulation is started with zero velocity, a uniform density n H,0 = 5 cm −3 , a uniform magnetic field B 0 = 4.6 µG and a uniform temperature T 0 = 10 4 K. The first few SN explosions rapidly bring the mean thermal, magnetic and kinetic energy to approximately steadystate values, with the magnetic field amplified to an rms value of 7.2 µG. The value of the mean magnetic field is chosen to achieve near equipartition with the kinetic energy at large scales, as shown in §4. It also yields an average value of |B| of 6.0 µG, consistent with the value of 6.0 ± 1.8 µG derived from the 'Millennium Arecibo 21cm Absorption-Line Survey' by Heiles & Troland (2005).
The simulation is run for 56 Myr (with a total of 359 SN explosions), initially without tracer particles and selfgravity. At t = 33 Myr we include 150 million passivelyadvected tracer particles following the mass distribution in the computational volume (each particle represents a fluid element of 0.013 M ). Tracer particles are advected with the same symplectic Kick-Drift-Kick scheme used for dark matter particles in Ramses, but the kicks are ignored -they are passive -and velocities in the drift are instead sampled using CIC interpolation of fluid velocities.
At t = 45 Myr we include self-gravity. Several clumps of order 1,000 M start to collapse. Larger resolution and sink particles are required at that point, which will be the topic of a followup work. For the purpose of this work, we stop the simulation after approximately 11 Myr of evolution with self-gravity. This is also motivated by the need to trace the exact location of SN explosions when clouds have been forming massive stars for a time of order 10 Myr, as explained in the next section.
To enforce a sub-pc resolution of the evolution of SN bubbles and of their interaction with the ISM, we adopt AMR criteria based on density, density gradients and pressure gradients. Although our root grid contains only 128 3 cells, our AMR criteria result in rather large volume filling factors of high-resolution cells: 75% of our computational domain is covered at a resolution equivalent to 256 3 , 22% at a resolution of 512 3 , and 2% at 1024 3 . Because the volume filling factor of the clouds selected in this work is approximately 0.5%, our clouds are all resolved at the maximum resolution.
The left panel of Figure 2 shows the logarithm of the projected density of the whole computational volume, at the final time of the simulation. The structure is highly filamentary and appears to be self-similar. The right panel of Figure 2 shows a sub-region magnified by a factor of four; the same type of filamentary structure is seen at this smaller scale. This structure is very similar to that previously found in supersonic simulations of isothermal turbulence, and consistent with the appearance of nearby GMCs mapped with the Herschel satellite.
Apart from SN explosions, the simulation neglects any other energy source, such as winds and radiation feedback from massive stars, or external forces, such as spiral arm shocks and fluctuations of the large-scale gravitational potential, that may affect MC turbulence in real galaxies. This choice allows us to test if SN turbulence alone can explain the origin and maintenance of GMC turbulence, while also providing a significant reduction of the computational cost. We also neglect the chemical evolution of the ISM. Although we do not expect Logarithm of projected density along the x-axis of the simulation volume. The mean magnetic field direction is horizontal on the image. The column density value is larger in darker regions, with black set to a maximum value of 5 × 10 22 cm −2 , and white to a minimum value of 6 × 10 20 cm −2 . The actual values of column density span a much larger range, and these limits have been chosen to optimize the contrast of the density structure. The image includes the whole computational volume, that is a size of 250 pc. Right Panel: Same as in the left panel, but for a region four times smaller (62.5 pc) shown by the white squared in the left panel. The color-table limit values have been increased by a factor of 4 and 5, 2 × 10 23 cm −2 (black) and 3 × 10 21 cm −2 (white), to match the larger mean density in this region.
that the dynamics would be strongly affected by a more precise computation of cooling and heating based on dynamically evolved chemical abundances, the selection of MCs and their comparison with observations would certainly benefit from a dynamical computation of the H2 and CO abundances (Glover & Clark 2012;Walch et al. 2014). The neglect of chemistry in this work is solely motivated by considerations of computational cost.

SN Explosions
Individual SN explosions are implemented with an instantaneous addition of 10 51 erg of thermal energy and 15 M of gas, distributed according to an exponential profile on a spherical region of radius r SN = 3dx = 0.73 pc, where dx is our smallest cell size, dx = L box /1024 = 0.24 pc. Kim & Ostriker (2015) have derived a useful condition for numerical convergence of the evolution of SN remnants, which states that both the grid size and the initial SN radius must be smaller than one third of the shell-formation radius, dx, r SN < r sf /3. Because in our case r SN > dx, and given the expression for r sf in Kim & Ostriker (2015), the condition for our simulation is r SN < 10 pc n −0.46 0 , where n 0 is the ambient density. Given our value of r SN = 0.73 pc, the condition becomes n 0 300 cm −3 . The volume filling factor of gas at density above that value varies in the approximate range between 0.0003 and 0.0009, with the largest value reached only at the end of the simulation, while the total number of SN explosions in the simulation is 359. Because the locations of SN explosions are randomly distributed, there is only a small probability that one or more SN explosions violate the condition for numerical convergence.
SN explosions are generated at random positions and times, neglecting the possibility of spatial and temporal clustering. In the galactic fountain simulations by de Avillez & Breitschwerdt (2005) and Joung et al. (2009) it was assumed that 60% of the SN explosions occurred in clusters, and only the remaining 40% at random locations. Joung et al. (2009) assumed that the clustered SN population is distributed in clusters with up to 40 stars more massive than 8 M , N SN = 40, with a probability of cluster size following a power law, ∼ N −2 SN . All massive stars from a given cluster were assumed to explode at the same location, though at different times distributed in an interval of 40 Myr, the approximate lifetimes of the least massive stars to explode.
However, assuming that clusters have a typical velocity dispersion of the cold gas of order 10 km/s, and considering the average value of N SN = 10, the typical cluster would contribute on average 1 SN explosion every 4 Myr, with a separation of approximately 40 pc between consecutive SN explosions, due to the cluster motion not accounted for by Joung et al. (2009). The 10 explosions would cover a distance of 400 pc over 40 Myr. Thus, even assuming the reasonable cluster statistics of Joung et al. (2009), the explosions would not appear to be clustered. A realistic prediction of spatial and temporal clustering of SN explosions requires that the formation and kinematics of individual stars is numerically resolved through the adoption of very high spatial resolution and sink particles. In the absence of that, a random SN distribution is a reasonable assumption.
Another potentially important issue is the correlation between the location of SN explosions and the position of their parent clouds. Recent simulations have shown that the large-scale structure of the ISM is sensitive to the amount of correlation between SN positions and density peaks Walch et al. 2014). Walch et al. (2014) have concluded that the random distribution case or the case of clustering in space and time, but not correlated with density peaks, are favored by the observations. On the contrary, in Dobbs & Pringle (2013) and Dobbs (2015) all SN explosions are assumed to occur inside MCs. As an MC (or any density peak above a certain threshold) is created, SN feedback is turned on in its interior. This unrealistic SN feedback does not allow to address the question of the origin of MC turbulence, as this is driven directly by SNe by design. Iffrig & Hennebelle (2015) have addressed this issue by focusing on the effect of a single SN on a single MC of 10 4 M , showing that the effect of a SN explosion on a MC is much stronger if the explosion occurs inside the cloud than outside of it. Their results may be affected by their specific choice of the ambient density, 1.2, 20 and 700 cm −3 , for the three SN positions they tested. The probability of high density is certainly increased inside a MC than outside of it. However, the probability of a SN explosion at high density must be quite small, because the volume filling fraction of dense filaments and clumps is low even within the volume enclosing a MC. Furthermore, HII regions and winds would probably prevent SN explosions in dense gas in general. Nevertheless, it seems reasonable that SN remnants expanding from within a MC would be more effective at bringing material above the escape velocity than SN remnants pushing on a MC from the outside.
With a Salpeter IMF (Salpeter 1955), the mean stellar mass between 8 and 100 M is approximately 19 M , corresponding to a stellar lifetime of approximately 9 Myr. Assuming that it takes at least 1 Myr to initiate star formation in a young cloud, and to accrete the stellar mass, the SN feedback would thus be important in a MC after a characteristic time of 10 Myr (assuming that massive stars can remain in the general region of their parent cloud for a parent-cloud crossing time). We will show in §8 that our most massive clouds of order 10 5 M selected towards the end of the simulation were formed approximately 20 Myr earlier (t form ≈ 20 Myr in Figure 21). For these clouds, thus, the SN feedback from locally-formed massive stars should start to play an important role, depending on the delay between cloud formation and the formation of the first massive stars.
The same is true for all clouds with lifetime of order of 10 Myr or larger. Because we find in §8 that the cloud lifetime is on average four times longer than the cloud dynamical time, and using the expression (25) for the dynamical time derived in §10.4, the condition is that the dynamical time is at least 5 Myr or longer, or, equivalently, that the cloud mass is larger than 500 M . As the formation of massive stars most likely requires clouds more massive than 500 M as well, we conclude that SN feedback from locally-formed massive stars should generally play a role for most clouds forming massive stars. As mentioned in §2, this is in fact the main reason why the simulation was stopped approximately 11 Myr after introducing self-gravity, as neglecting the locally-formed massive stars beyond that time would be unrealistic. The effect of the correlation of SN explosions with their parent clouds will be considered in a future work where the for-mation and kinematics of individual massive stars will be numerically resolved in order to model self-consistently the precise time and location of SN explosions.

MC SELECTION
The main question addressed by this work, whether SN explosions can drive and maintain the observed MC turbulence, can be partly answered independently of the definition of MCs, by computing the velocity structure functions for the dense and cold gas throughout the computational volume. As long as most of such gas is in clouds, its global structure functions should be equivalent to the average of the structure functions of MCs. Nevertheless, it is important to characterize the turbulence within individual clouds and to also compute other cloud properties that can be compared with observations.
Because chemistry is not included in this work, we can only define MCs as cold over-densities in the SN-driven ISM turbulence. A detailed comparison with the observations is beyond the scope of this work. It would require synthetic observations, hence molecular abundances from chemical-network calculations. In the absence of synthetic observations, we prefer to avoid the selection of MCs in position-position-velocity (PPV) space, and instead define MCs in 3D space (PPP) in the simplest possible way, as connected regions above a single threshold gas density, n H,min . In order to test if our results depend on spatial resolution or threshold density, we extract clouds using three different mesh sizes, 2dx, 4dx and 8dx (we create uniform grids of 512 3 , 256 3 and 128 3 cells, respectively), and four different threshold densities, n H,min = 100, 200, 400 and 800 cm −3 , generating 12 cloud catalogs. Examples of clouds from one such catalog (n H,min = 100 cm −3 and 512 3 resolution) are shown in Figure 3. The images show the projected density of 18 clouds, the 2nd to the 7th most massive ones in each of three snapshots.
The analysis of MC properties derived with the full three-dimensional information, such as the results on MC velocity scaling, the discussion on the virial parameter and cloud structure, the evaluation of cloud lifetimes and the study of the cloud magnetic field, is based on clouds extracted with n H,min = 100 cm −3 and 128 3 resolution. For this analysis, the resolution only affects the definition of cloud boundaries, as all results are derived from the position, velocity, density and magnetic field of the tracer particles, thus taking advantage of the highest resolution.
When we compare with observational data, deriving projected quantities such as surface density, equivalent radius and line-of-sight velocity dispersion, we verify the results on all 12 catalogs. For the mass and size distributions, we also report quantitatively on the dependence of the slope of their power-law tails on cloud extraction density and resolution. However, all the plots we show in the comparison with the observations are based on the highest-resolution catalog (512 3 cells, or 0.49 pc) and on the threshold density that best matches the observed mass-size relation, n H,min = 200 cm −3 . Furthermore, velocity dispersions are based on tracer particles and so take advantage of the highest available resolution, dx = 0.24 pc. 2 This resolution matches well the high-  , integrating over the whole simulation volume, independent of density. The horizontal dotted line shows the initial magnetic energy (also the magnetic energy corresponding to the mean magnetic field in the simulation, which is a conserved quantity). There is near equipartition between the lowest values of kinetic energy, the mean value of the magnetic energy, and the highest values of the thermal energy.
est one in the observational survey we consider (Heyer et al. 2001). After selecting clouds with circular velocity v c < 20 km s −1 and mass M cl > 100 M , we are left with 3,228 observed clouds with measured distances corresponding to a range of spatial resolutions of 0.24 − 3.0 pc.
We select from the simulations only clouds more massive than 100 M to guarantee that, even at the lowest value of n H,min = 100 cm −3 , the smallest clouds contain more than 1,000 computational cells (our largest cloud of nearly 3 × 10 5 M contains more than 3×10 6 cells). Because we use 150 million tracer particles, our smallest and lowest-threshold-density clouds contain a minimum of approximately 7,000 particles, while our most massive cloud of 3 × 10 5 M contains more than 20 million particles. With this minimum mass, each simulation snapshot yields over 200 clouds. In the comparison with the observations (see §10) we use 7 snapshots from approximately the last 6 Myr of the simulation, to include the effect of self-gravity, resulting in sample sizes ranging from 595 clouds in the smallest catalog (128 3 resolution and n H,min = 800 cm −3 ), to 1,615 in the largest one (512 3 resolution and n H,min = 100 cm −3 ). The plots we show in §10, based on 512 3 resolution and n H,min = 200 cm −3 , use measurements from 1,547 clouds. The catalog sizes could be considered three times larger, as projected quantities are computed in the three orthogonal directions. However, the clouds (and cloud masses) would be the same, so the three samples would not be completely independent. Thus, we compare with the observations using only one of the directions, after verifying that the results are independent of direction.
Despite the mass range of over three orders of magnitude of our clouds, we refer to all of them as MCs, instead of following the common nomenclature that classifies clouds as cores, clumps, molecular clouds, giant molecular clouds, and giant molecular cloud complexes,  Figure 4, but with the energies computed only for gas with density above 100 cm −3 , the same threshold density adopted for the cloud selection. There is a clear energy separation in the dense gas, with the lowest kinetic energy values being approximately an order of magnitude larger than the mean magnetic energy, and two orders of magnitude larger than the lowest values of the thermal energy. Thus, the turbulence in the dense gas is both supersonic and super-Alfvénic.
in order of increasing mass. As we view all clouds as cold density enhancements of the ISM turbulence, and because of the scale-free nature of the turbulence, that nomenclature is not useful for this work.

TOTAL ISM ENERGIES
We analyze our simulation in a time interval starting after the turbulence has been fully developed until the end of the simulation, between t = 33.06 Myr and t = 56.43 Myr. This is also the time interval for which we have tracer particles in the simulation. The evolution of the total kinetic, magnetic and thermal energies in that time interval is shown in Figure 4. The three energies are all around 10 51 erg. Both the thermal (E th ) and kinetic (E k ) energies show strong oscillations, while the magnetic energy (E m ) is almost constant with time. The kinetic energy oscillates mostly above 10 51 erg, and, interestingly, at its valleys, the magnetic energy is almost in equipartition with E k . The thermal energy is the smallest, mostly below both the magnetic and kinetic energies.
Given the strong temporal oscillations in kinetic energy, it may seem surprising that the magnetic energy does not experience significant fluctuations at all. The reason is that the kinetic energy peaks correspond to the energy injected by SN explosions mainly in the form of expansions and shocks. While the magnetic field can be strongly compressed and amplified by the shocks, the volume filling factor of the postshock gas is tiny, and thus the total magnetic energy of the computational volume is barely changed, as illustrated by the following estimate.
Consider a SN remnant of radius R SN and a preshock magnetic field B 1 . Compression by the SN shock would amplify the magnetic field strength by a factor equal to the density-jump factor, Γ, of the shock. In the adiabatic phase, Γ = 4 and in the radiative phase, Γ M A with M A = v sh /(B 1 / √ 4πρ 1 ) the Alfvénic Mach number of the shock. As the remnant expands, solenoidal turbulent motions also develop (see §5.1), which can amplify the magnetic field as well. The timescale for the magnetic energy amplification by solenoidal motions is roughly the turnover time of the largest eddies in the postshock region (e.g. Federrath et al. 2011a). We assume that the turbulent velocity in the postshock region is the postshock velocity, v sh /Γ, and the large eddy size is the thickness of the postshock region of the remnant, which is R SN /(3Γ) (due to the compression by the shock). The large eddy turnover time is then estimated to be R SN /(3v sh ). Considering the age of the remnant is 2 5 R SN /v sh and 2 7 R SN /v sh in the adiabatic and radiative phases, respectively, solenoidal motions may amplify the magnetic energy by one e-fold or so, meaning an amplification factor of A 3 for the magnetic energy. Including both the effects of shocks and solenoidal motions, the magnetic energy density in the postshock region would be 1 8π AΓ 2 B 2 1 . The above description of the magnetic field amplification only applies to the compressed layer behind the shock. Due to the compression by a factor of Γ, the width of the compressed layer is given by R SN /(3Γ), so the total magnetic energy within the compressed layer is estimated to be 1 6 AΓB 2 1 R 3 SN . On the other hand, the magnetic energy in the hot cavity interior to the compressed layer is small due to the expansion and can be neglected. Considering that the magnetic strength is still B 1 in regions not reached by the SN remnant, the total magnetic energy in the simulation box is where N SN is the number of SNe exploded around the same time and V box is the volume of the simulation box. Because the highest kinetic-energy peaks are of the order 10 52 erg, N SN may be as large as 10. If we define a filling factor of each SN remnant as f = 4πR 3 SN /(3V box ), the total magnetic energy can be written as 1 8π [N SN (AΓ − 1)f + 1]B 2 1 V box , meaning that E m is amplified by a factor of N SN (AΓ − 1)f + 1. In the Sedov phase, Γ is constant, and the amplification factor for the magnetic energy increases with the filling factor, f . Applying the physical conditions in our simulation, we find that, when the Sedov phase ends due to radiative cooling, f increases to 3 × 10 −4 . Thus, with Γ = 4 and A 3, (AΓ − 1)f is only 0.003. Therefore, due to the tiny filling factor of the SN remnant, the total magnetic energy is amplified only by a negligible amount during its early evolution.
When the remnant evolution enters the radiative phase, the jumping factor, Γ M A , which can be significantly larger than 4. Using the pressure-driven snow-plough solution to account for the deceleration of the shock velocity (and hence the decrease of M A with time) and the increase of R SN (and the filling factor, f ), we find that the maximum of the amplification factor, (A 2 Γ − 1)f , by a single SN in the radiative phase is 0.05. Therefore, even if at a given time there are 10 SNRs that happen to simultaneously amplify the magnetic energy by the maximum possible amount, the total amplification of the magnetic energy is only 50%. Note that, in the above estimate, we have ignored the dynamical effect of the magnetic field and the dissipation of the magnetic energy, and thus the realistic amplification due to the SN explosions may be considerably smaller than 50%. This explains the near constancy of the total magnetic energy despite the strong oscillations in the total kinetic energy of the flow. Although the SN energy is introduced purely as thermal energy in the simulation, the thermal energy is efficiently converted into kinetic energy, and is also partly radiated, so the resulting turbulence is mildly supersonic, with a time-averaged kinetic to thermal energy ratio of E k /E th = 3.75. A crucial question for this work is what fraction of that total kinetic energy is given to the dense gas and, thus, is available as turbulent kinetic energy for the MC turbulence. The answer is given in Figure 5, showing the time evolution of the total energies in the gas with density larger than 100 cm −3 . The time-averaged ratio of dense-gas kinetic energy and total kinetic energy is E k,d /E k = 0.11. The ratio of the kinetic energies per unit mass is (E k,d /M d )/(E k /M box ) = 0.55 (the densegas mass fraction is M d /M box = 0.18), and oscillates in time between 0.19 and 0.97, with the largest values achieved during the peaks of total kinetic energy. This large ratio means that, per unit mass, the kinetic energy is only a factor of ∼ 2 less than equally distributed between the dense gas and the rest of the gas. This high efficiency of kinetic energy transfer to the dense gas results in a realistic velocity dispersion of dense gas, σ v,d = 8.5 km/s, as shown in Figure 6. Figure 5 also shows a clear separation of energies in the dense gas, with E k,d /E m,d = 27.4 and E m,d /E th,d = 9.7. Thus, the turbulence in the dense gas is both supersonic and super-Alfvénic, as further confirmed below for individual clouds selected from our simulation (see §10) and first suggested by , 1999.

POWER SPECTRA AND DRIVING SCALE
The expansion of SN remnants deposit energy on a broad range of scales, affecting the velocity scaling of the turbulence. We observe this in the form of wave-like features in the velocity power spectrum or deviations in the velocity structure functions (see §6) at small and intermediate scales in nearly one fourth of the snapshots of our simulation, as expected for our rate of approximately 6 SN explosions per Myr, and a characteristic expansion time to 20-40 pc of ∼ 4 × 10 4 yr. Figure 7 shows the power spectra from a single snapshot that captures the early expansion phase of a single SN remnant (here it has reached a diameter of approx-  (k), square root of kinetic energy, E K (k), and magnetic field, Em(k), obtained from the average of the power spectra of the x, y and z components of these fields in the root grid (128 3 computational cells) of a single snapshot at time t = 52.75 Myr. This time captures the early expansion of a SN remnant to a diameter of approximately 30 pc, as indicated by the peak of the velocity power spectrum. The wavy appearance of Ec(k) is explained in the text. The fact that Es does not show similar fluctuations, and Es < Ec at the energy peak, indicates that this remnant has so far expanded into a relatively uniform hot medium, where the baroclinic effect is negligible and solenoidal modes are not efficiently generated yet.
imately 30 pc) that has not had any major interaction with dense gas yet, so most of the freshly-injected energy is still in the compressive modes (see §5.1) 3 . The Figure shows the power of the velocity, E(k), of the solenoidal (divergence-free) velocity component from a Helmholtz decomposition, E s (k), of the compressive (curl-free) velocity component from the same Helmholtz decomposition, E c (k), of the square root of kinetic energy (ρ 1/2 u i ), E K (k), and of the magnetic field, E m (k). They are computed from the root grid (128 3 computational cells) of a single snapshot at time t = 52.75 Myr. The velocity power spectra, especially the compressive one, exhibit significant fluctuations at large and intermediate k. The wavy behavior is a signature of strong SN shocks and can be understood with a one-dimensional illustration. Consider two SN shocks moving in opposite directions away from the origin (the explosion center). If the velocity profile in between the two shocks is assumed to be roughly linear with the distance from the origin, it is straightforward to show that the power spectrum is ∝ k −2 (cos(kR) − sin(kR)/kR) 2 , where R is the radius of the SN shock. Clearly, this spectrum oscillates at k R −1 , which explains the wavy behavior of the velocity spectra shown in Figure 7.
In between SN explosions, or in regions not directly affected by a rapidly expanding SN remnant, the flow should have time to relax from the transient state with direct SN impact, and we expect the velocity spectra to be smoother. To test this, we consider a time range t = 40.5−48.8 Myr, around the middle of our integration time. This choice avoids the larger SN rate towards the  Figure 7, but averaged over 66 snapshots within the time interval t = 40.5−48.8 Myr. During this time interval the random SN rate is a bit lower than during the first and last third of the simulation, so there is a reduced probability that a snapshot captures the early expansion of SN remnants with the corresponding perturbations to the power spectra. As a result, the time average is representative of the statistically relaxed power spectra. The power-law slopes are obtained from a least-square fit in the range of wave numbers 4.2 ≤ k L box /2π ≤ 12.1.
beginning and the end of the simulation (see Figure 4), and thus minimizes the number of snapshots with direct impact from very recent SN explosions. Figure 8 shows the average three-dimensional power spectra computed from the root grid of 66 snapshots in the chosen time range. For each snapshot, the 3D spectra are obtained by averaging the three components (e.g., the average of the power spectra of u x , u y and u z ). The average spectra we obtain are qualitatively similar to those of the most relaxed snapshots, and they exhibit a clear (though short) inertial range. The slopes have been computed in the wavenumber interval 4.2 ≤ k L box /2π ≤ 12.1, corresponding to the scale interval 59.3pc ≥ ≥ 20.6pc, where all spectra are very well described by power laws. The measured inertial-range slopes are given in Figure 8. To our knowledge, this is the first time that inertial-range slopes are identified in SN-driven turbulence. Because of the importance of this result, a full discussion of the power spectra will be presented elsewhere. Here, we focus only on two points of direct interest for this work, the ratio of compressive to solenoidal modes and the driving scale.

Ratio of Compressive to Solenoidal Modes
SN driving should not be viewed as a purely compressive form of driving. Our simulation indicates that, around the effective driving scale (see §5.2), compressive modes are not dominant over solenoidal ones. We find that the compressive-to-solenoidal ratio is typically smaller than unity, as shown by the time-averaged power spectra in Figure 8. 4 After a SN explosion, the momentum of the ejecta is gradually transferred to the ISM during the expansion of the SN remnant until the final momentum-conserving snow-plough phase. As a result, SN-driving covers a wide The time t 0 is just before the SN explosion, while the times t 1 and t 2 are after the explosions, when the diameter of the remnant is approximately 40 and 80 pc, respectively. Between the times t 1 and t 2 , the expansion of the SN remnant brings the power in both the compressive and solenoidal modes to larger scales, as explained in the text. Right: Squared root of projected temperature (upper row of panels) and logarithm of projected density (lower row of panels) for the three times at which the power spectra are computed. At time t 1 the edge of the remnant is barely visible in the logarithm of projected density, while at time t 2 the gas density within the remnant has significantly decreased in the lowest density regions.
range of scales and persists for a relatively long time after the explosion. This complex driving process cannot be primarily compressive because vorticity is readily generated by the baroclinic effect and amplified by the non-linear advection, as explained below. Even the ISM forcing immediately after the explosion is far from purely compressive, as hydrodynamical instabilities of the blast wave start already in the interior of the star, generating vorticity and making the ejecta very clumpy and asymmetric (e.g. Chevalier & Klein 1978;Herant & Woosley 1994;Kifonidis et al. 2006;Couch et al. 2009;Wongwathanarat et al. 2015). In our simulation, the forcing at the moment of the thermal energy injection is not purely compressive either, because the acceleration from the corresponding pressure force has a non-zero curl (a baroclinic term) due to the non-uniformity of the medium (see the argument below).
To examine the ratio of compressive and solenoidal modes, it is helpful to write down the equations of the flow divergence and vorticity, and, The terms on the right-hand sides arise from the pressure term in the Navier-Stokes equation, and (∇ρ × ∇p)/ρ 2 in Equation (2) is called the baroclinic term, which contributes to generate vortical motions when the gradients of ρ and p are not aligned. The (ω · ∇)u term is known as vortex stretching and can amplify the vorticity. In our simulation, the flow velocity is driven by the pressure term. If we denote as p s the contribution to the pressure from the SN thermal energy source, the effective driving acceleration is −(∇p s )/ρ. It follows immediately from Equations (1) and (2) that the divergence and curl of the effective acceleration are given by (∇ρ · ∇p s )/ρ 2 − (∇ 2 p s )/ρ and (∇ρ × ∇p s )/ρ 2 , respec-tively. Because the SN locations are selected randomly in our simulation and due to the density fluctuations in the ambient gas, the pressure and density gradients (∇ρ and ∇p s ) around the boundary of the initial SN sphere are not aligned in general, meaning that, at the instant of the SN energy injection, the effective driving acceleration at the boundary of the SN sphere is not curl-free. Therefore, the effective driving in our simulation is not purely compressive.
If (∇ρ · ∇p s )/ρ 2 dominates over (∇ 2 p s )/ρ (which is supported by the fact that |∇ ln(ρ)| · |∇p| > |∇ 2 p| in a significant fraction of our simulated flow), the ratio of compressive to solenoidal power generated by the pressure source is determined mainly by the angle, θ, between the directions of ∇ρ and ∇p s . Because in our simulation the SN locations are random, one may expect that the angle between the gradients is also random, so that (cos θ) 2 1/3 and (sin θ) 2 2/3. In that case, the effective driving generates the vorticity variance (corresponding to solenoidal motions) twice faster than the divergence variance. Although of a qualitative nature, the above argument provides an explanation of why the solenoidal spectrum is typically larger than the compressive one.
By monitoring solenoidal and compressive spectra as a function of time, we can observe that the expansions of the SN remnants bring both the compressive and solenoidal modes to larger scales, as shown in Figure  9. The power in compressive modes moves to larger scales according to the simplified 1D model given earlier, which also explains the wavy features in the spectrum; the solenoidal power is transferred to larger scales by the expansion through the non-linear term in Equation (2), ω(∇ · u). Figure 9 shows the velocity power spectra in a (100 pc) 3 volume centered around a SN explosion. The spectra are computed at three different times, one just before the SN explosion, and two after the SN remnant has reached an approximate size of 40 and 80 pc. The velocity field is multiplied by a tapered-cosine window function before performing the Helmholtz decomposition in Fourier space, to gradually set the velocity to zero at the volume boundaries, making the velocity field periodic. The transfer of both solenoidal and compressive power to larger scale is clearly seen between the times t 1 and t 2 , as the remnant expands from 40 to 80 pc.
At times or in regions not too close to very young SN remnants, the solenoidal and compressive modes have a chance to develop cascades toward small scales and to interact with each other, evolving toward a dynamically and statistically relaxed state. The timeaveraged spectra in Figure 8 indicate that in the relaxed state the solenoidal component of the velocity is larger than the compressive component at all k. The ratio, At the smallest wave numbers, the modes are almost in equipartition, χ(k) ∼ 1. However, E s (k) is remarkably shallow, while E c (k) ∝ k −2 , so the compressiveto-solenoidal ratio decreases rapidly towards larger wave numbers, χ(k) ∝ k −0.67 in the inertial range. Although they did not obtain clear power-law scaling and did not identify inertial-range slopes in their SN-driven simulation, Balsara et al. (2004) also found that χ(k) decreases with increasing k and χ(k) 1 at large wave numbers. The different inertial-range slopes of E c (k) and E s (k) and the rapid decrease of χ(k) with increasing k are in contrast to previous studies, which typically found that χ(k) is more or less constant in the inertial range. For example, simulations of weakly compressible turbulence with Mach numbers M s ∼ 1 showed that both compressive and solenoidal modes have Kolmogorov-like inertial-range spectrum, E c (k) ∝ E s (k) ∝ k −5/3 , and χ(k) ∼ 0.05 (e.g. Porter et al. 1998Porter et al. , 1999Porter et al. , 2002. At the opposite extreme, simulations with purely compressive driving find χ(k) ∼ 1 in the inertial range, and the same Burgers-like slopes for both modes, E c (k) ∼ E s (k) ∝ k −2 (Federrath 2013). Simulations with solenoidal (or mixed) driving and large Mach number, M s 1, yield values χ(k) ≈ 0.3 − 0.5 Federrath 2013), with only a slight decline towards large wave numbers.
An important difference between these studies and our simulation is that they all adopted a barotropic equation of state, assuming either adiabaticity or isothermality. The baroclinic effect is thus absent in all the simulations mentioned above. As discussed earlier, when SNe explode in our simulation, the baroclinic effect from the pressure source immediately drive solenoidal motions around the SN spheres at a rate similar to compressible modes. As the SN remnant and the flow evolve, the general baroclinic effect in Equation (2) is also likely to play an important role as the flow develops a cascade and evolves toward relaxation. As compressive motions in the form of shocks or expansions encounter a dense region in the flow, the baroclinic effect gives rise to shear and vortices around and within the region. Since the baroclinic term depends on the gradients of the density and pressure, the conversion from compressive to solenoidal motions is expected to be more efficient at smaller scales. This explains why the compressive-to-solenoidal ratio decreases rapidly with increasing k in our simulation. We thus argue that the existence of the baroclinic effect in our simulation is responsible for the different behaviors of E c (k), E s (k) and χ(k) compared to previous barotropic simulations.
There has been recent interest in the regime of supersonic turbulence with purely compressive driving, showing that it affects the probability distribution of gas density Molina et al. 2012), the magnetic field amplification by turbulence (Federrath et al. 2011a), and the inertial-range velocity scaling Federrath 2013), with possible consequences for models of star formation or for turbulence in galaxy clusters (Porter et al. 2015). However, our result that SN-driving is not primarily compressive, particularly at MC scales, suggests that features specific to isothermal flows with highly compressive driving may not apply to ISM turbulence. If SN shock-waves are the main energy source for MC turbulence, the driving acceleration would consist of a significant fraction of solenoidal modes, which arise through the baroclinic effect, when the SN shock impacts a cloud. Thus, idealized isothermal simulations with purely solenoidal driving may better capture MC turbulence than isothermal simulations with purely compressive driving, and previous results on turbulent fragmentation based on solenoidal driving may not require significant correction. A careful study of the full implication for star formation of our result would require the evaluation of the compressive-to-solenoidal ratio specifically for the clouds formed in the SN-driven turbulence, which we pursue in a separate work (Pan et al. 2015).

Energy Injection Scale
Using the velocity power spectrum, E(k), we can define a length, L in , that corresponds approximately to the scale where most of kinetic energy is contained, and can thus be interpreted as a characteristic scale of energy injection by SN explosions: The time dependence of L in is shown in Figure 10, where we have also plotted the rms velocity, σ v . The value of L in oscillates between approximately 50 and 100 pc, with many of the peaks corresponding also to peaks in σ v . The time average and standard deviations are 70.5 pc and 12.0 pc respectively. Figure 10 also shows the transverse integral scale, L 22 . The longitudinal and transverse integral scales, L 11 and L 22 , are defined as the integrals of the two-point velocity correlation functions: where R L and R N are the longitudinal and transverse correlation functions, with u and u n being the velocity component parallel to and one (of the two) transverse component perpendicular to , respectively. L 11 and L 22 are typically smaller than the injection length, L in . In isotropic turbulence, exact relations exist between L 11 , L 22 and L in . For example, in incompressible turbulence, L 11 = 2L 22 = 3L in /8 (Monin & Iaglom 1975). Our simulated flow is highly compressible, and we find that the longitudinal and transverse integral scales are close to each other with L 11 = 19.9 pc, and L 22 = 19.4 pc. The near equality of L 11 and L 22 suggests that kinetic energies contained in solenoidal and compressive modes are comparable. Figure 8 shows that at the scale of L in the solenoidal and compressive spectra are in equipartition, in the sense that the solenoidal spectrum is twice larger than the compressive one due to the extra degree of freedom. One can demonstrate that in such case L 11 and L 22 should be exactly equal, as found in our simulation. Furthermore, in that case L 11 and L 22 are expected to be equal to L in /4, which further explains the ratios, L in /L 11 = 3.5 and L in /L 22 = 3.6. By comparing four galactic-fountain simulations with different SN rates, Joung et al. (2009) found that the energy injection scale decreases with increasing SN rate. In the model with approximately the same SN rate as in our simulation they obtained L in = 87 pc. This value is approximately consistent with the one derived here, if we account for the fact that Joung et al. (2009) assumed that 60% of the SN explosions are spatially correlated, as discussed above in §2.1, which enhances the formation of super-bubbles and so should tend to increase the correlation length.
Using their galactic-fountain simulation, de Avillez & Breitschwerdt (2007) computed the longitudinal and transverse integral scales, L 11 and L 22 . They found that L 22 /L 11 = 0.5 − 0.6, consistent with the ratio in incompressible turbulence, implying that the overall compressibility of their simulated flow was likely low. Their measured value of 75.2 pc for L 11 would indicate an injection length scale of L in 200 pc. This injection length scale is significantly larger than L in 70 pc found in our simulation. Based on the dependence of the integral scale on the SN rate from Joung et al. (2009), this difference may be attributed to the much smaller SN rate in the simulations analyzed by de Avillez & Breitschwerdt (2007).

VELOCITY STRUCTURE
In order to characterize the turbulence, we compute the velocity structure functions, first for the whole computational volume, and then for individual clouds. The velocity structure functions of order p are defined as: where the velocity component u is parallel (longitudinal structure function) or perpendicular (transversal structure function) to the vector and the spatial average is over all positions x. Boldyrev (2002) proposed an extension to supersonic turbulence of the intermittency model by She and Lévêque (She & Lévêque 1994;Dubrulle 1994): This velocity scaling has been found to provide a very accurate prediction for numerical simulations of highly supersonic and super-Alfvénic turbulence Padoan et al. 2004b;Pan & Scannapieco 2011). Padoan et al. (2004b) showed that, as the rms Mach number of the turbulence increases, the structure function scaling varies from the She-Lévêque scaling of incompressible turbulence to the Boldyrev's scaling, which has been interpreted as a gradual change of the Hausdorff dimension of the most dissipative structures from 1 (dissipation in filaments) to 2 (dissipation in sheets). Although the scaling of equation (7) may not apply to flows driven by purely compressive forces ), we have shown in §5.1 that the compressive modes are not dominant in SN-driven turbulence.
Because of the limited extent of the turbulence inertial range in numerical simulations, the structure functions are usually power laws only over a very limited range of scales, if at all. Thus, the scaling exponents are usually derived by normalizing the structure functions to the third-order one, which yields power laws that extend well into the (numerical or physical) dissipation range of scales. This useful property is known as 'extended self-similarity' (Benzi et al. 1993).
The third-order structure function always yields a slope larger than unity in supersonic turbulence, while ζ(3) = 1 is an exact result in incompressible turbulence, the so called '4/5 law' first derived by Kolmogorov (1941). This is because in supersonic turbulence the third-order structure functions should be computed with some density-weighting factor. For example, a density weight inspired by the assumption of constant energy transfer of a compressible flow, u ∼ ( /ρ) 1/3 (Fleck 1996), was proposed by Kritsuk et al. (2007) and further tested by Kowal & Lazarian (2007) and by Schmidt et al. (2008). In this work, we compute the structure functions either using the tracer-particle positions and velocities, or using the gas velocity values on a uniform mesh, with a density-weighting method that is equivalent to the use of tracer particles. This density-weighting method is different from that proposed by Kritsuk et al. (2007).
To calculate the structure functions based on tracer particles in a MC, we search particle pairs at given dis- Fig. 11.-Density-weighted longitudinal velocity structure functions, S dw p ( ) (see equation (9)), of orders p = 1, 2 and 3, plotted versus the separation , computed over the whole computational volume (filled circles). The plots show the time average of the structure functions from 24 snapshots covering uniformly the time interval t = 33 − 56 Myr (one snapshot per Myr). The solid lines corresponds to the structure function slopes predicted by Boldyrev (2002) and confirmed in simulations of randomly-driven isothermal turbulence Padoan et al. 2004b;Pan & Scannapieco 2011), ζ(1) = 0.42, ζ(2) = 0.74, ζ(3) = 1.0. Although they provide an excellent fit, these solid lines are not obtained by fitting the data.
tances, , compute their relative velocities, and average over all the particle pairs to obtain the relative velocity moments at different orders, p. For example, the p-th order structure function is computed as where N is the number of pairs separated by a distance of , and u 1n − u 2n is the relative velocity of the n-th pair of particles. These tracer-based structure functions have a built-in density weighting, because the number of particle pairs at a given distance depends on the flow densities at the particle positions (the tracer particles are initialized with a number density proportional to the local gas density). For structure functions over the entire simulation box, we compute grid-based structure functions with a density weighting defined as It is straightforward to prove that S tr p ( ) and S dw p ( ) are statistically equivalent to each other, if the number of tracers is large enough to sufficiently sample the flow density field. To show this, consider two infinitesimal volumes, dV 1 and dV 2 , around two points, 1 and 2, separated by a distance of . When computing the tracerbased structure functions, these two volumes would contribute to a total number of N 12 = (n p /ρ) 2 ρ 1 ρ 2 dV 1 dV 2 particle pairs at a distance of , wheren p andρ are the mean particle number density and mean flow density, respectively. In other words, the two points 1 and 2 are essentially counted N 12 ∝ ρ 1 ρ 2 times in the computation of S tr p ( ). This is equivalent to a density-weighting factor of ∝ ρ 1 ρ 2 , suggesting that S tr p ( ) is equal to the grid-based structure function, S dw p ( ). The density-weighting scheme adopted here is of particular interest in the light of a result derived by Falkovich et al. (2010). Motivated by Kolmogorov's 4/5 law for incompressible flows, Falkovich et al. (2010) obtained an exact relation for compressible turbulence, where p is the pressure of the flow. In highly supersonic turbulence, the pressure term may be neglected, and we have The quantity ρ(x)ρ(x + )u i (x)u i (x + )u j (x + ) is closely related (although not exactly equivalent) to the density-weighted third order structure function, S dw 3 . We therefore expect a linear scaling for our density-weighted third order structure function, S dw 3 ∝ , which turns out to be confirmed by our simulation for both the entire flow and individual MCs.

Global Velocity Scaling
We first compute the velocity structure functions averaged over the whole computational volume. The velocity field is first remapped onto a uniform mesh of 512 3 cells, and then the structure functions are computed with the density-weighting method described above. The volume contains large voids of very low density, with a very low number of tracer particles, so we prefer this mesh-based method instead of the equivalent tracer-particle method (which we use below for the structure functions inside MCs). To speed up the calculations, we only consider velocity differences along the three main orthogonal directions, and 16 values of cell distances, . Even with this limitation, the sample size is large enough to yield reliable low-order statistics. Figure 11 shows the first, second, and third order longitudinal velocity structure functions, S dw p ( ), plotted versus the separation , computed over the whole volume as described above, and time-averaged using 24 snapshots covering uniformly the time interval t = 33 − 56 Myr (one snapshot per Myr). They are well approximated by power laws in the approximate range of scales To compare with values previously found in studies of randomly-driven supersonic isothermal turbulence, in this and following figures we over-plot lines showing the slopes from eq. 7. Notice that eq. 7 only gives the exponents of the velocity structure functions without density weighting and normalized to the third-order one. Thus, in this comparison, we make the reasonable assumption that the exponents normalized to the third order should be the same for the unweighted structure functions as for the density-weighted ones. However, future works should recompute the velocity structure functions of randomlydriven, supersonic isothermal turbulence simulations using the same density-weighting scheme as in this work.
While the first and second order structure functions in Figure 11 have exponents within 3 sigma and 1 sigma respectively from the values of equation 7, ζ(3) is significantly larger than unity, despite the time averaging. We find that this deviation from unity, and even from a power law, occurs in snapshots during periods with the highest SN rate, when it is more likely that a snapshot is very close in time to a very recent SN explosion. During the brief, initial period of the SN bubble expansion, SN driving has a direct effect on a range of scales. On the contrary, in snapshots that are not too close in time to SN explosions, the SN remnants have already expanded to large scale, the turbulence has had time to relax, and the third order structure function is found to be a nearly perfect power law with ζ(3) = 1. Figure 12 shows an example of velocity structure functions from a single snapshot, at t = 44.5 Myr, that is not too close in time to a SN explosion. The third order structure function is now a power law, with the expected numerical decay below approximately 3-4 pc, and the third order slope is indistinguishable from unity. The exponents in this single snapshot are ζ(1) = 0.41 ± 0.01, ζ(2) = 0.77 ± 0.02, ζ(3) = 1.00 ± 0.02, consistent with equation 7 within one or two sigma. To better quantify the uncertainty of these exponents, and to illustrate their time dependence, we compute their values for each snapshot in the time interval t = 33 − 56 Myr (8 snapshots per Myr), by fitting the structure functions as a function of in the range between 3 and 30 pc, and plot them versus time in Figure 13. The two pairs of solid and dotted horizontal lines show the predicted values for incompressible (lower values) and supersonic (higher values) turbulence (She-Lévêque and Bodyrev scaling respectively). One can see that the values are usually within the predicted ones for the first and second order exponents, while the third order exponent is systematically above unity during the first and last thirds of that time interval, when the SN rate and the kinetic energy appear to be a bit higher than in the middle third (see Figure 4). The time average of the singlesnapshot values are ζ(1) = 0.39±0.04, ζ(2) = 0.75±0.05, ζ(3) = 1.1 ± 0.2, well within one sigma of the values from equation 7.
The logarithmic fluctuations in Figure 13 (the y-axis of the plot is logarithmic) are significantly larger for the third order than for the first and second orders, meaning that deviations from the expected scaling due to the effect of SN driving is stronger for higher order statistics. Therefore, it is not convenient to take advantage of extended self-similarity and measure the exponents by plotting the structure functions versus the third order one, instead of versus , at least not for the first and second order case, as their time evolution would then have much larger fluctuations, due to the larger fluctuations of the third order one. Joung & Mac Low (2006) and de Avillez & Breitschwerdt (2007) computed the velocity scaling exponents from galactic-fountain simulations and found a scaling law consistent with Boldyrev's prediction. However, because of insufficient dynamic range below the driving scale, their simulations did not yield a power-law inertial range, and so the exponents could be computed only relative to the third order one (taking advantage of the extended self-similarity). Therefore, their result cannot be compared directly with the ISM velocity scaling derived from observations. Furthermore, their simulations did not reach the necessary spatial resolution to describe the evolution of individual SN remnants and their interaction with MCs in detail, and to study the velocity scaling within individual clouds, to test if MC turbulence is consistent with SN driving.
Because our simulation yields power-law velocity structure functions as a function of , and the values of the exponents are the same as in previous numerical studies of MC turbulence based on random driving with a large-scale volume force, we conclude that the use of an artificial force in those previous studies did not result in incorrect velocity scaling, so no major corrections to IMF and SFR models based on turbulent fragmentation should be needed. However, because deviations from the average scaling laws are found in snapshots that are very close in time to SN explosions, it may be worthwhile investigating the possibility of minor effects on star formation resulting from such deviations.

Velocity Scaling within MCs
In order to test if SN driving can generate the same velocity scaling also inside MCs, despite their density contrast, we compute the velocity structure functions inside MCs selected from our simulation as described in §3. To better constrain the scaling exponents, we have computed the velocity structure functions for the 15 most massive clouds of each snapshot. Because of the very complex cloud shapes, we find it convenient to compute the structure functions based on the position and velocity of the tracer particles, S tr p ( ) (see eq. 8). MCs are regions of density enhancement, so their velocity field is sampled well by the tracer particles. In fact, they contain so many tracer particles that, to speed up the calculation, we randomly select a number of particle pairs 500 times smaller than all possible pairs in each cloud, resulting in approximately 2 to 200 million pairs per cloud. As a further simplification, we also compute the differences of each velocity component irrespective of their orientation relative to the separation vector, , so we do not distinguish between transversal and longitudinal structure functions.
We find velocity differences growing with scale up to 10-30 pc, depending on the cloud size. As an example, Figure 14 shows the average of the structure functions of the 15 most massive clouds selected at time t = 34.2 Myr. The structure functions inside these GMCs are consistent with the global ones and are well approximated by power laws down to a scale of approximately 10dx = 2.4 pc. At a separation = 4dx = 0.96 pc, the velocity differences are only approximately 30% below the value extrapolated from the inertial-range scaling, as shown by the first order structure function in Figure 14.
In Figures 15 and 16 we show examples of structure functions of individual MCs, one with no evidence of the direct effect of SN driving (Figure 15), and one with significant deviations at ≤ 1 pc, due to a recent nearby SN explosion ( Figure 16).
These results show that, despite the large contrast between the MC density and the average ambient density, the kinetic energy of SN explosions is effectively transferred into turbulence within individual clouds, where it establishes the usual scaling laws of supersonic turbulence, all the way to the smallest scales where the simulation is affected by numerical dissipation. Real MCs also exhibit power-law velocity scaling (e.g. Heyer & Brunt 2004;Padoan et al. 2006). The slope of ζ(2) = 0.8 ± 0.1 derived for the Perseus region by Padoan et al. (2006) using the method by Lazarian & Pogosyan (2000) is formally consistent with the scaling laws derived here for the MCs of our simulation. On the other hand, the principlecomponent-analysis results by Heyer & Brunt (2004) give a very large slope, ζ(2) = 1.12 ± 0.04, which is hard to interpret, as it is steeper than the scaling from the Burgers equation that models an infinitely compressible flow. A careful study of the consistency between the structure function slopes of the SN-driven simulation and of the observations requires the computation of synthetic observations and is beyond the scope of this work.
All the plots in this section adopt an arbitrary normalization of the structure functions. The actual normalization of MC turbulence, that is the velocity dispersion of clouds of a given size, is discussed below, comparing the velocity-size Larson relation of our clouds with that from the observations. Larson (1981) interpreted the MC velocity-size relation he discovered as due to a turbulent cascade in the ISM. Because of the very large Reynolds number of the observed motions in the cold ISM, the development of a turbulent cascade is unavoidable, and both analytical arguments and numerical simulations have demonstrated that Larson relations can be viewed as the natural consequence of supersonic turbulence (e.g. Kritsuk et al. 2013). Nevertheless, the velocity-size relation has also been interpreted as the consequence of the MC self-gravity (e.g. Solomon et al. 1987;Heyer et al. 2009), because MC virial masses are often comparable to the masses estimated from the CO luminosity.

VIRIAL PARAMETER AND CLOUD STRUCTURE
The relative importance of turbulence and self-gravity is measured by the virial parameter, introduced by Bertoldi & McKee (1992): where σ v is the one-dimensional velocity dispersion, and the dynamical time is defined as: The last equality in (12) is exact in the case of an idealized spherical cloud of uniform density. For more realistic cloud mass distributions, the virial parameter is only an approximation of the ratio of kinetic and gravitational energies.
To estimate the relative importance of turbulence and self-gravity in clouds from our simulation, we have computed the virial parameter and the kinetic and gravitational energies of clouds selected from six snapshots, three before and three after the inclusion of gravity in the simulation, using a threshold density n H,min = 100 cm −3 . The virial parameter of a cloud is measured using the positions and velocities of the tracer particles belonging to that cloud, and defining the cloud radius as the  Figure 17, showing that the average ratio of α vir and 2E k /Eg and its scatter are nearly unchanged for most clouds, relative to their value before the inclusion of self-gravity. Clouds with very low values of 2E k /Eg contain collapsing cores whose gravitational energy has been included in the computation of the cloud Eg. They are star-forming clouds with relatively low value of α vir , but do not show evidence of global collapse.
rms of the particle positions: wherex i ≡ N n=1 x i,n /N are the components of the mean particle position (the cloud's barycenter) and N is the total number of tracer particles in the cloud.
While the virial parameter only depends on global MC properties (mass, radius and rms velocity), the ratio of kinetic to gravitational energy is sensitive to the cloud internal structure (mass distribution and shape) and to correlations between density and velocity (e.g. Federrath & Klessen 2012). Thus, the comparison between α vir and 2E k /E g is a useful tool to probe the internal structure of MCs and its evolution under the effect of self-gravity. As in the case of the virial parameter, we compute E k and E g of a cloud using the velocities and positions of the tracer particles in that cloud: where N is the total number of tracer particles in the cloud, m is the mass associated to a tracer particle, u i is the modulus of the velocity of the i-th particle, and r ij is the distance between the i-th and the j-th particles. This expression for E g assumes that the cloud is isolated, as in the definition of the virial parameter. In future work, this expression should be contrasted with a formulation that accounts self-consistently for the surrounding mass distribution by using the actual gravitational potential from the simulation (e.g. Federrath & Klessen 2012). -Probability distributions of the viral parameter for the same three snapshots before self-gravity as in Figure 17 (shaded histogram) and for the three snapshots after self-gravity as in Figure  18 (unshaded histogram). Figure 17 shows the comparison between α vir and 2E k /E g for clouds selected from three snapshots covering a time span of 4 Myr before the inclusion of selfgravity. Figure 18 shows the same plot, but based on three snapshots after the inclusion of self-gravity, covering the last 4 Myr of the simulation. On average, the clouds of Figure 18 are selected 9 Myr after the inclusion of self-gravity, while their average free-fall time, based on the density estimated as n H,cl = M cl /(4πR 3 cl /3) (see Figure 20), is 4.4 ± 2.1 Myr, and their average dynamical time, t dyn ≡ R cl /σ v,3D , is 2.6 ± 1.5 Myr. Thus, self-gravity has been active for approximately two cloud free-fall times and four cloud dynamical times on average, with significant changes occurring only in regions that have too small filling factors (e.g. small collapsing cores) to change global statistics significantly. Indeed, gravity is not expected to be important for clouds with α vir > 1. Selecting clouds with α vir ≤ 1.0, 0.5 and 0.25, we find a mean free-fall time of 2.1 ± 1.3 Myr, 1.5 ± 1.3 Myr and 1.2 ± 1.0 Myr respectively. Thus, for the clouds where it should be most important, self-gravity has been active for 4 to 8 free-fall times on average. To establish whether the absence of significant trends during this time interval will continue for the lifetime of the various structures requires, as discussed in Section 2.1, future simulations that include the supernova feedback from the massive stars produced by the simulation itself. Figure 17 shows that α vir provides a remarkably good approximation of the energy ratio, 2E k /E g , despite the complexity of the cloud structure. We find the ratio α vir /(2E k /E g ) = 1.20 on average, constant over three orders of magnitude in E k /E g (we have verified that it is constant also over the full range of cloud masses) and with a small scatter (the standard deviation is approximately 20% of the mean); it is also nearly unchanged after gravity is included in the simulation. Figure 18 shows that most clouds follow approximately the same relation of α vir versus 2E k /E g as in the case without self-gravity, except that a few of them have significantly lower values of 2E k /E g , because they contain collapsed cores that have also been included in the computation of the total E g . The comparison of the two figures, without and with gravity, shows the main effect of self-gravity  Figure 17 (shaded histogram) and for the three snapshots after self-gravity as in Figure  18 (unshaded histogram).
is to cause the collapse of dense cores inside the clouds (hence star formation), while the cloud virial parameter is not strongly affected.
In Figure 19 we show the probability distribution of α vir for the three snapshots without gravity (shaded histogram) and with gravity (unshaded histogram), where we have included also the star-forming clouds with very low values of 2E k /E g . The two distributions are very similar, with the one including gravity slightly shifted towards lower values; the mean values are 8.5 without gravity and 6.6 with gravity. The small shift between the two distributions shows that self-gravity causes some amount of cloud contraction, but not a significant change in global cloud structure. This is further confirmed by the histograms of cloud density shown in Figure 20, where the cloud density is defined as n H,cl = M cl /(4πR 3 cl /3). The histogram for the clouds with self-gravity is only slightly shifted to higher density, with the mean value changing from 183 cm −3 to 264 cm −3 , before and after the introduction of gravity respectively. This small increase in cloud density shows that self-gravity does not cause a global cloud collapse, even if it drives star formation through the collapse of dense cores within MCs. Only clouds with α vir 10 contribute to star formation, according to Figure 18, with most star formation occurring in clouds with α vir 3, in agreement with recent studies of the SFR in supersonic turbulence, showing that the SFR is mainly controlled by the virial parameter Padoan et al. 2012). Future simulations, where self-gravity is active for much longer than the 11 Myr period of our run, and where selfconsistent supernova feedback is included, will be needed to further test the above results.

MC LIFETIMES
While the virial parameter estimates the relative importance of turbulence and self-gravity, the comparison of the cloud lifetime with either the dynamical time or the free-fall time provides a definitive assessment of the actual dynamical state of clouds. For example, shortlived clouds are evidently not gravitationally bound, while their virial parameter may be of order unity and thus would not allow to draw a definitive conclusion about their dynamical state.
Being so difficult to constrain observationally, MC lifetimes are a valuable outcome of numerical modeling. They can be measured in our simulation thanks to the introduction of tracer particles, but only up to a maximum age of 23 Myr, the time interval with tracer particles in our simulation (11 Myr if we considered only clouds during the time interval with self-gravity). Although this maximum age is comparable to or smaller than the lifetime of the largest clouds, it is at least significantly longer than the mean free-fall time (4.4 Myr) and the mean dynamical time (2.6 Myr) of the clouds in the simulation, so it allows us to evaluate the influence of gravity on the cloud lifetimes for most clouds in our study. Dobbs and Pringle (2013) measured MC lifetimes by defining the continuation of a cloud at a later time as that with the largest number of particles in common with that cloud. The time when the number of particles (and the mass) in common drops to half or less of that in the original cloud marks the end of the cloud lifetime. Cloud precursors and formation times are defined in the same way, by considering the clouds at earlier times with the largest number of particles in common with the original cloud. As pointed out by the authors, this method has the drawback that clouds precursors or continuations are often not found with more than half of the original mass because of changes in density contours, rather than a real cloud dispersion. For example, based on this method, a cloud may seem to have dispersed after a certain time, but may later reappear.
After verifying in our own simulation that this lifetime definition is indeed very uncertain, we have chosen to estimate cloud lifetimes with a different method that is independent of the specific density contours of clouds in past and future snapshots. Because the formation of a cloud implies converging flows (even in the absence of self-gravity), and its dispersion requires diverging flows, we simply define the lifetime of a cloud as the time interval during which the cloud radius, defined always by all the tracer particles belonging to that cloud, is within a factor of two from the radius at the time the cloud is selected. In other words, we follow the cloud tracer particles in the past, until their radius has doubled in size, which marks the formation time of the cloud, and in the future also until the radius has doubled in size, which marks the dispersion time of the cloud. Figure 21 shows the formation time, t form (interval between cloud formation and time of cloud selection), and dispersion time, t disp (interval between time of cloud selection and cloud dispersion) versus the cloud mass, using all the clouds more massive than approximately 700 M from each of 12 snapshots covering the whole time interval with tracer particles (the time between snapshots is 2 Myr). The sum of the two times gives the cloud lifetime, t life = t form + t disp , but we have plotted the two times separately (empty squares and circles) because in many cases we can identify the cloud formation time and not its dispersion time (for clouds selected towards the end of the simulation), or vice-versa (for clouds selected shortly after the introduction of tracer particles). The average values of t form and t disp are biased by the large number of upper limits (not plotted to avoid confusion). In order to obtain a nearly unbiased estimate of the average values, we consider only clouds selected from the first and the last snapshots of the series (filled circles and filled squares, respectively, in Figure 21). The first snapshot has 12 clouds more massive than 700 M , yielding 10 measured values of t disp and only two upper limits; the last snapshots has 40 clouds more massive than 700 M , 39 of which with a measured value of t form and only one with an upper limit to t form . Based on these measurements alone, we find t disp = 11.1 Myr and t form = 9.2 Myr. By assuming that clouds are selected at a random moment of their lifetime, we should expect t life = 2 t form = 2 t disp , and so we can average together all the 49 values and obtain t life = 21.4 Myr.
This lifetime is based on clouds covering over two orders of magnitude in mass, and only on 39 measurements. We can further improve our estimate of MC lifetimes by using all measurements and, at the same time, derive the mass dependence of the lifetimes, if we normalize the lifetime to the cloud dynamical time, t dyn . Figure  22 plots t life versus t dyn of 64 clouds for which we have measured values of both t form and t disp (filled circles), besides the formation and dispersion times for all the other clouds where these are measured (empty squares and circles respectively). The plot in Figure 22 shows an approximately linear correlation between t life and 4 t dyn . The mean and standard deviation of the ratio of the two times are t life /t dyn = 4.5 ± 1.4.
As shown by Figure 23, this estimate is derived from clouds with M cl 6 × 10 3 M , and thus it should be considered only as an extrapolation when applied to more massive MCs. Nevertheless, this result for the t life /t dyn ratio is consistent with the measured formation times of the most massive MCs in the simulation. Figure 23 shows that the four most massive GMCs, with masses around 10 5 M , have t form /t dyn close to two, which would imply t life /t dyn ≈ 4, in the absence of selection biases due to the limited time interval covered by the tracer particles. Furthermore, while we cannot use all the values of t form and t disp from Figure 21 to estimate an unbiased average lifetime from t life ≈ 2 t form ≈ 2 t disp , we can still use all of the corresponding ratios, t form /t dyn and t disp /t dyn , to estimate an unbiased average ratio of lifetime to dynamical time from t life /t dyn ≈ 2 t form /t dyn ≈ 2 t disp /t dyn , if we assume that this ratio is independent of cloud mass. Indeed, the dashed and dotted lines in Figure 23 show that t form /t dyn ≈ t disp /t dyn ≈ 2, using values over the whole mass range, consistent with the estimate t life /t dyn ≈ 4, based on clouds with M cl 6 × 10 3 M . This result shows that both the formation and the dispersion of the MCs in our sample take two dynamical times, on average. This is an indication that both the formation and the dispersion of the MCs in our sample is controlled by the turbulence, with little influence of self-gravity. Because of the non-negligible scatter in the ratio of cloud lifetime to dynamical time, one may expect that at least the clouds with the largest ratios may have longer lifetime due to their self-gravity. This is not the case: we have verified that there is actually a positive correlation between t life /t dyn and α vir ,meaning that larger values of t life /t dyn are usually due to smaller values of t dyn because of larger σ v (hence larger α vir ), rather than longer t life as a consequence of a lower α vir . Thus, there is no significant imprint of self-gravity in the cloud lifetimes, even if more than half of our clouds are selected at a time after self-gravity has been included in the simulation. Future simulations, where selfconsistent supernova feedback allows longer runs with selfgravity, are needed to test if this lack of significant imprint of selfgravity continues, and if it extends to MCs with longer lifetimes and to higher surface density MCs.
We should also stress the caveat that, for the most massive MCs of ∼ 10 5 M we could only measure formation times, and not dispersion times, due to their long lifetime (and dynamical time) and the limited duration of the simulation. Thus, we cannot rule out that, at least the most massive MCs, could have dispersion times significantly longer than two dynamical times. However, that would imply dispersion times longer than 20 Myr (life-  Figure 22. The long-dashed line shows the mean ratio for the clouds with both dispersion and formation times measured (filled circles), (t form + t disp )/t dyn = 4.5. The dotted and the dashed lines the mean ratio for clouds with only formation or dispersion times respectively, t form /t dyn = 2.2, t disp /t dyn = 2.4. All these average values are consistent with t life /t dyn ≈ 4, independent of cloud mass.
times longer than 40 Myr) for such clouds, a timescale over which the extra energy injection from SN explosions of locally formed massive stars would presumably succeed in dispersing the clouds, even if the general ISM turbulence could not (see discussion in §2.1).
To derive actual values of cloud lifetime as a function of cloud mass, taking advantage of our result (17)

MAGNETIC FIELD IN MCS AND MC FORMATION
Our simulation adopts a mean magnetic-field strength consistent with the Galactic one (see §2), so the magnetic field inside clouds selected from the simulation should be comparable to that in real MCs. We have already shown in Figures 4 and 5 that the mean magnetic energy is not far from equipartition with the mean thermal and kinetic energies averaged over the whole volume, while the energy ratios are much larger in the dense gas. This clear energy separation in dense gas, with E k,d /E m,d = 25.1 and E m,d /E th,d = 9.8, is the necessary consequence of the near equipartition at the largest scales. Being only mildly super-Alfvénic, large-scale compressive motions cannot compress the mean magnetic field by a large factor, so the density enhancement of MCs is largely achieved with compressions along field lines, resulting in a mean magnetic field strength in the dense gas not much larger than the total mean field. The mean magnetic field of 4.6 µG is amplified by the SN-driven turbulence to an rms value of 7.2 µG, averaged over the whole volume and between t = 33 and 56 Myr. The rms field strength in the dense gas is 12.8 µG, not even a factor of two larger.
To investigate the role of the magnetic field in individual MCs, we consider the same catalog of 1,547 clouds as in the comparison with the observations discussed in the next section. The clouds are selected from 7 snapshots during the final 6 Myr of the simulation, at a resolution of 0.49 pc and with a density threshold of n H,min = 200  Figures 31 and 34. The dashed line shows the mean magnetic field averaged over the whole computation volume (also the initial mean field). Empty circles correspond to the mean value of the magnetic field of all tracer particles in each cloud, while filled circles give the rms value. cm −3 . We compute both the mean and the rms magnetic field of each cloud using the values sampled by the tracer particles, B = Σ i B i /N and B 2 1/2 = (Σ i B 2 i /N ) 1/2 , where B i is the magnetic field strength sampled by the particle i in a given cloud, and N is the total number of particles in that cloud. These magnetic field values are plotted versus cloud mass in Figure 24, where the horizontal dashed line represents the mean magnetic field in the computational volume, B 0 = 4.6 µG. The mean field in the clouds is approximately 10 µG on the average, only twice larger than B 0 , and independent of cloud mass. We have verified that the mean magnetic field strength of the clouds is also independent of their mean gas density.
The relatively small increase of the cloud mean magnetic field relative to B 0 and its independence of gas density are characteristic of trans-Alfvénic supersonic turbulence , 1999, and further illustrates that MCs must be formed by compressive motions primarily along magnetic field lines, due to the non-negligible magnetic pressure prior to the compression and cooling of the low-density gas. As the gas is being compressed into a nascent MC by random large-scale motions, the increasing density and decreasing cooling time cause a drop in both the Alfvén and sound speeds . As a result, the turbulence within a MC is super-Alfvénic and highly supersonic, while the larger-scale flows responsible for its formation are trans-Alfvénic and mildly supersonic. Because in super-Alfvénic turbulence the magnetic field is amplified by compressions, as shown by a positive B − n correlation , 1999, dense cores formed by shocks within MCs (Padoan et al. 2001) have an enhanced magnetic-field strength on average. Furthermore, cores are topologically the ultimate zero-dimensional destination of a fluid element undergoing compression, as they can be viewed as the intersection of filaments that are formed by the intersection of postshock sheets. Much of the flow turbulent energy is dissipated by the time it 'stagnates' into a core. Due to this drop in turbulent energy, together with the increase in magnetic-field strength and density, the turbulence inside dense cores  Figure 24, computed with the cloud mean magnetic field (empty circles), or the cloud rms magnetic field (filled cirlces).
is trans-Alfvénic and trans-sonic. In summary, super-Alfvénic and supersonic MC turbulence is the natural consequence of large-scale trans-Alfvénic trans-sonic turbulence and also the natural origin of small-scale trans-Alfvénic trans-sonic turbulence in prestellar cores.
The small-scale enhancement of the magnetic field within MCs is partly illustrated in Figure 24 by the values of the rms field in the clouds that is approximately a factor of two larger than the mean field in the most massive clouds. As a more direct demonstration of the super-Alfvénic nature of MC turbulence, Figure 25 shows the cloud rms Alfvénic Mach number versus the cloud mass. The Mach number is computed as the ratio of the cloud rms velocity and the cloud Alfvén velocity, where the latter is computed either with the mean magnetic field (empty circles) or with the rms magnetic field (filled circles), and using the mean density sampled by the tracer particles. Nearly all clouds with mass larger than 10 3 M are super-Alfvénic, even considering their amplified field strength. For the 41 GMCs with masses larger than 10 4 M , the average Alfvénic Mach number is 8.3 with respect to the mean field, and 3.9 with respect to the rms field.
The super-Alfvénic nature of the turbulence in the clouds from our simulation is consistent with the observational evidence. Based on the comparison between simulations of MHD turbulence and MC observations, , 1999 suggested that MC turbulence was better characterized by supersonic turbulent flows with M A 1 than flows with M A ≈ 1. This result was later confirmed with the aid of synthetic observations (Padoan et al. 2004a) and synthetic Zeeman splitting measurements (Lunttila et al. 2008(Lunttila et al. , 2009. Taking advantage of the anisotropy of MHD turbulence, Heyer & Brunt (2012) demonstrated that the densest regions of the Taurus MC complex are characterized by super-Alfvénic turbulence, while in low density regions the motions are sub or trans-Alfvénic, also consistent with the picture from our simulation, where MCs are formed by large-scale trans-Alfvénic turbulence, and thus fed preferentially by motions along magnetic field lines, as discussed above (Nordlund & Padoan 2003;Padoan et al. 2010). Because the clouds are selected as regions with density n H > 100 cm −3 , no tracer particles belonging to a cloud has density lower than that. To extend the relation to densities n H < 100 cm −3 , we define a cloud volume delimited by the smallest and largest coordinates of the tracer particles of that cloud, and include the values of B and n H of each cell of that volume to compute the relation. The two solid lines are least-square fits for densities n H < 10 3 cm −3 and n H > 10 3 cm −3 , giving B ∼ n 0.13 H and B ∼ n 0.29 H , respectively.
To further characterize the cloud turbulence, we have computed the B − n relation inside all the clouds of our sample, using again the values of B and n of the tracer particles inside the clouds. We divide the density in logarithmic bins, and compute the mean magnetic field strength and its standard deviation in each density bin. We then average these values among all the clouds, using weights proportional to the number of tracer particles in the density bins of each cloud. Figure 26 shows this mean B − n relation. We also illustrate the large scatter by plotting error bars that correspond to twice the standard deviation above and below the mean values. The two solid lines are least-square fits for n H < 10 3 cm −3 , B ∼ n 0.13 H , and for n H > 10 3 cm −3 , B ∼ n 0.29 H . The stronger dependence of the magnetic-field strength on density at n H > 10 3 cm −3 than at lower density is qualitatively consistent with the observations (Crutcher et al. 2010). The slope we derive is much smaller than that derived by Crutcher et al. (2010) at high densities, B ∼ n 0.65 H . However, their slope does not refer to the mean magnetic field at a given density, but to its maximum value. Considering the large number of measured upper limits well below such maximum values, the dependence of the mean B on density could be significantly shallower than the estimated slope of the upper envelope of the B − n relation. Furthermore, the Bayesian analysis by Crutcher et al. (2010) assumes a uniform distribution of the magnetic field strength, while this distribution is exponential in super-Alfvénic turbulence (Padoan & Nordlund 1999) (we have verified it is exponential also in our clouds). Finally, and most importantly, the B − n relation in Figure 26 is not computed for a selection of dense cores, as in the observations, but using every single tracer particle in the cloud, so it should not be compared quantitatively to the observational B − n relation. Such a comparison would require synthetic Zeeman observations of a selection of dense cores, as in Lunttila et al. (2008Lunttila et al. ( , 2009. It would also require higher spatial resolution, because most of the observed cores with the largest detected magnetic-field strengths, at densities of order 10 5 − 10 7 cm −3 , have sizes substantially smaller than the spatial resolution of our simulation. The higher resolution would also allow to better resolve the dynamo amplification in dense cores (Federrath et al. 2011b), which would tend to increase the slope of the B − n relation.

COMPARISON WITH OUTER-GALAXY MCS
To further test our results, we carry out a comparison of the properties of our MCs with those of observed MCs. This is a preliminary approach based on the derivation of projected quantities, such as column density, equivalent radius and line-of-sight velocity dispersion. Follow-up studies with synthetic observations taking into consideration chemistry and radiative transfer are also needed.
Our observational sample of choice for this comparison is the MC catalog by Heyer et al. (2001), extracted from a decomposition of the 12 CO FCRAO Outer Galaxy Survey . Besides the large dynamic range of the survey, its main advantage is that for the Outer Galaxy there is no blending of emission from separate MCs along the line of sight, or at least the problem is strongly mitigated compared with the Inner Galaxy. As a result, a very large number of clouds can be reliably selected over a broad range of cloud masses and sizes. The catalog contains a total number of 10,156 objects, up to a mass of approximately 8 × 10 5 M and a size of 45 pc. It is estimated to be complete down to a mass of approximately 600 M and a cloud size of 3 pc.
Inner-Galaxy MC catalogs are far less reliable and complete because of velocity blending, so they are not suitable for the comparison we pursue here. For example, the recent catalog of Inner-Galaxy MCs Roman-Duval et al. 2009) extracted from the UB-FCRAO Galactic Ring Survey (Jackson et al. 2006) contains objects between 1 and 10 6 M , but is estimated to be complete only above 4 × 10 4 M (Roman-Duval et al. 2009). We suspect a more realistic completeness limit may be 2 × 10 5 M , because the mass distribution is a power law only above that mass (Roman-Duval et al. 2009), which corresponds to a size limit of approximately 20 pc, judging from the mass-size relation and from the lack of a power law in the size distribution below 20 pc. This severe incompleteness suggests that the cloud surface density may be overestimated by a large factor. The MC mass distribution is expected to be a power law down to small masses, so the abrupt departure from a power law below 2 × 10 5 M indicates that much of the missing mass from smaller clouds is incorrectly assigned to larger ones due to velocity blending. The failure to select individual three-dimensional clouds is also demonstrated by the absence of a velocity-size correlation and, possibly, by the extremely low values of the cloud virial parameters (the distribution peaks at α vir ≈ 0.2), which would imply a larger star formation rate and a stronger signature of global collapse than observed. In summary, the differences between Galactic-Ring and Outer-Galaxy MCs may not be as large as often assumed. Thus, although our comparison is primarily with Outer-Galaxy MCs, the  Heyer et al. (2001). The dotted-line histogram is the mass distribution of the observational sample, excluding clouds with radius Re ≤ 3.1 pc, the size-completeness limit of the Outer-Galaxy survey. It shows a contribution to the mass distribution by clouds below this completeness limit up to a mass of approximately 1,000 M . The leastsquare-fit to the observational mass distribution above 1.5 × 10 3 M has a slope of −0.99 ± 0.09 and is shown by the dashed-dotted line.
results may be applicable to Galactic clouds in general.
As in Heyer et al. (2001), we consider only the subset of 3,901 clouds with circular velocities v c < −20 km s −1 , because of kinematic distance accuracy. We further select clouds with mass M cl > 100 M , as that is the mass limit for our numerical cloud catalogs, resulting in a total sample of 3,228 Outer-Galaxy MCs. Given the distances to the clouds and the angular resolution of the survey, the spatial resolution varies between 0.4 pc and 3.8 pc. Therefore, the cloud extraction of our highest resolution catalog with dx = 0.48 pc matches well the highest resolution in the observations. The main limitation of the survey is the velocity resolution, only slightly better than 1 km s −1 , which, combined with the measurement of line width based on the equivalent width instead of the antenna-temperature-averaged velocity dispersion, results in a minimum velocity dispersion of clouds of approximately 0.5 km s −1 . However, we show that the data can be used to test both the slope and the normalization of the Larson velocity-size relation from the simulation despite this low velocity resolution.
A explained in §3, we illustrate this comparison using clouds selected from 7 snapshots from the last 6 Myr of the simulation. However, we have verified that all the observational MC properties discussed in this section are essentially the same when derived from a catalog of clouds selected from the last 6 Myr prior to the inclusion of gravity. This confirms that global MC properties are primarily the result of SN-driven turbulence, with little modification due to self-gravity, apart from the slight increase in mean cloud density shown in §7, and an increase in the mass of the largest cloud.
Of the 12 cloud catalogs described in §3, we choose the one with the highest-resolution (512 3 cells, or 0.49 pc) and with the threshold density that best matches the observed mass-size relation, n H,min = 200 cm −3 , for all the plots in this section, and we compute velocity dispersions using the tracer particles, thus taking advantage of the highest resolution of the simulation, dx = 0.24 pc, due to the large number density of tracer particles within dense clouds. This catalog contains 1,547 objects. Figure 27 shows the mass distribution of our clouds (shaded histogram). The histogram is well approximated by a power law, with a slope of −0.88 ± 0.06 in the approximate mass range 200 − 10 5 M (dashed line). The figure also shows the mass distribution of the MCs from the observational sample that is also nicely fit by a power law, with a slope of −0.91±0.09 in the approximate mass range 1.5 × 10 3 − 2 × 10 5 M (dashed-dotted line). Heyer et al. (2001) derived a slope of −0.80 ± 0.03 by including all clouds down to the completeness limit of 600 M . Our slope is a bit steeper because we are a bit more conservative on the completeness limit. The dotted-line histogram in Figure 27 shows the mass distribution for a sub-sample where we include only clouds above the size completeness limit of 3.1 pc, the value derived by Heyer et al. (2001). The comparison with the histogram of the full sample shows that some MCs with sizes below the size completeness limit are found with masses up to approximately 1,000 M , so we consider the catalog to be complete only above that mass value.

Mass Distribution
The mass distribution of our clouds is consistent with that of the MCs from the Outer Galaxy Survey. Furthermore, this is true for the mass distribution from all numerical cloud catalogs we have compiled, independent of the numerical resolution and threshold density of cloud selection. As shown in Figure 28, the exponent of the power-law fit to the mass distribution has a weak dependence on resolution and threshold density and, within the 1-σ error bars shown in the plot, it is consistent with the observational exponent in all cases.
The largest cloud in our n H,min = 200 cm −3 catalog Fig. 29.-Probability distribution of cloud size for the same cloud catalog from the simulation (shaded histogram) and the same observational sample (unshaded, solid-line histogram) as in Figure  27. The dotted-line histogram is the size distribution of the observational sample, excluding clouds with mass M cl < 600 M , the mass-completeness limit of the Outer-Galaxy survey. The dashed line is a fit to the tail of the histogram from the simulation, with slope −2.5 ± 0.2, and the dashed-dotted line a fit to the observational size distribution, with slope −2.3 ± 0.3. has a mass of 1.3 × 10 5 M , a few times smaller than the largest MC in the observational sample. However, Figure 27 shows that this maximum mass is consistent with the largest mass expected from our sample size and the slope of the mass distribution. If we simulated a region larger than 250 pc, and thus collected a much larger cloud sample, we would likely derive a power-law mass distribution extended to larger masses. The observations are consistent with a power-law mass distribution for clouds more massive than those in the Outer Galaxy Survey. For example, in their analysis of the MC sample by Solomon et al. (1987), Williams & McKee (1997) found a comparable slope of −0.81 ± 0.14 in the mass range 3 × 10 5 − 5.6 × 10 6 M . Although more uncertain, the slope they obtained from the analysis of the sample by Scoville et al. (1987), −0.67 ± 0.25, is also consistent with the slope of the Outer-Galaxy MCs at lower masses. More recently, Roman-Duval et al. (2010) found a slope of −0.64 ± 0.25 in the mass range 4 × 10 4 -10 6 M from the analysis of the Galactic Ring Survey (a more conservative completeness limit of 8 × 10 4 M gives a slope of −0.86 ± 0.25). Following Heyer et al. (2001) and most observational works, as a measure of a cloud's size we adopt the equivalent radius, R e ≡ (A cl /π) 0.5 , where A cl is the cloud projected area. The probability distribution of R e for the clouds from the simulation is shown in Figure 29 (shaded histogram). It is well approximated by a power law with a slope of −2.5 ± 0.2 in the approximate range of 2 − 15 pc. Within the uncertainty, this is consistent with the slope of −2.3 ± 0.3 of the observational sample in the approximate equivalent-radius range of 5 − 50 pc. Furthermore, Figure 28 shows that the size distributions of clouds selected from the simulation with different threshold density and resolution are also consistent with the observations, within the 1-σ uncertainty (except for the catalog with n H,min = 200 cm −3 and the lowest-resolution). As in the case of the mass distribution, we have been slightly more conservative in the estimation of the size completeness limit, based on evidence that around the value of 3.1 pc, the completeness limit estimated by Heyer et al. (2001), we still find some contribution from clouds with mass below the mass completeness limit of 600 M , as illustrated by the dotted histogram in Figure 29. As a result, we find the same slope as in Heyer et al. (2001), but with a three times larger uncertainty. The same power law seems to apply to even larger clouds. For example, Sanders et al. (1985) find a slope of −2.3 ± 0.2 for a sample of 80 clouds in the approximate size range of 20 − 80 pc.

Size Distribution
Because we have previously computed a three dimensional cloud radius, R cl , from the simulation, we can test the relation between the observable radius, R e , and the three-dimensional one. The comparison is shown in Figure 30, where the dashed line corresponds to R e = R cl . R e is smaller than R cl for most clouds with R cl 2 pc, and larger than R cl for most clouds with R cl 2 pc. The average ratio is R e /R cl = 0.87 and increases towards smaller radii. As a consequence, the probability distribution of R cl is slightly shallower than that of R e .

Velocity-Size Relation
The velocity-size relation of MCs determines the normalization of the velocity scaling inside individual clouds. In §6.2, we have shown that the energy from SN explosions sets a turbulent cascade inside individual MCs that follows the usual velocity scaling of supersonic turbulence, but we have not discussed the normalization of the velocity structure functions of individual clouds. To address the velocity normalization, we compute the internal rms velocity of our clouds based on the velocity of their tracer particles. Being derived from tracer particles, this rms velocity is mass-weighted, which is a reasonable approximation when comparing it with the antenna-temperature-weighted rms velocity from MC observations. As in the observations, we compute the onedimensional (line-of-sight) rms velocity, in the direction perpendicular to the plane (plane of the sky) where we measure the equivalent radius. Figure 31 shows the velocity-size relation of the clouds from the simulation (empty circles) and from the obser-  Figures 27 and 29, but the observations are shown only for Re > 4 pc, because for smaller cloud sizes the lower envelope of the velocity-size relation is not resolved by the observation (the minimum value of σv that can be detected is ∼ 0.5 km s −1 ). The thin solid line shows the mean values of σv in logarithmic bins of Re, and the thick solid line is a fit to those values, giving a slope of 0.39 ± 0.07. The dashed line is the fit to the binned data from the observations, giving a slope of 0.48 ± 0.06, and exactly the same normalization as the simulation at Re ≈ 10 pc.
vational sample (filled circles). We have excluded all the observed clouds with R e < 4 pc, because the observations cannot detect velocity dispersions smaller than approximately 0.5 km s −1 , due to the low velocity resolution. The lower envelope of the velocity-size relation of both the simulation and the observations decreases with decreasing cloud size, and reaches the value of 0.5 km s −1 at approximately 4 pc. Thus, the velocity dispersion of a fraction of the clouds smaller than 4 pc may be significantly overestimated, with that fraction growing towards smaller cloud radii, causing an artificial flattening of the velocity-size relation. As shown in Figure 6 of Heyer et al. (2001), the velocity-size relation for the full sample is essentially flat below 4 pc.
Values of σ v from the simulation are not to be trusted for the smallest cloud sizes, R e < 2 pc, because of the increasing effect of numerical dissipation towards smaller scales. The thin solid line in Figure 31 shows the average values of σ v in logarithmic bins of R e . While it is nicely fit by a power law for R e > 2 pc, it clearly drops at smaller cloud sizes. This is consistent with the cloud structure functions shown in Figure 14, where the numerical dissipation starts to become important below approximately 2 pc as well.
Despite these limitations imposed by the low velocity resolution of the observations and the numerical dissipation in the simulation, we still have a sufficient range in R e where the simulation and the observations can be compared. Both the upper and the lower envelopes of the velocity-size relation are very similar in the two cases. Furthermore, the velocity normalization is nearly identical. The thick solid and dashed lines in Figure 31 show the least square fits of the average values of σ v in logarithmic bins of R e of the simulation (for R e > 2 pc) and of the observations (for R e > 4 pc), respectively. From the simulation we get and from the observations: so the velocity normalization at R e = 10 pc is indistinguishable in the two cases. This agreement between the simulation and the observations in the slope, total scatter and normalization of the velocity-size relation is strong evidence that SN driving alone can be responsible for the turbulence observed in MCs. The universality of the MC velocity normalization has been questioned by Heyer et al. (2009), claiming that it depends on column density, and thus that the velocitysize relation is controlled by gravity rather than being a natural consequence of the ISM turbulence. To further confirm the agreement between the simulation and the observations, we show the velocity normalization as a function of column density in Figure 32 (we plot σ v /R 0.5 e as in Heyer et al. (2009), even if the slope of the velocity-size relation is actually smaller than 0.5). Because the range of cloud column densities is similar in the two cases, we should not expect a different normalization even if it depended on surface density. Figure 32 shows a good overlap between our clouds and the observations. We also plot the values for the Galactic-Ring clouds in Roman-Duval et al. (2010) that have an average surface density an order of magnitude larger than the Outer-Galaxy clouds (notice that the difference in surface density is only a factor of five for equal cloud mass or size, and could have been overestimated by a factor of two or three, as explained in the opening of §10). The figure shows that there is no difference in the velocity normalization of Galactic-Ring and Outer-Galaxy MCs, despite the difference in surface density. Thus, we conclude that the normalization of the velocity-size relation of the MCs in our sample is consistent with being controlled by SN-driven turbulence, rather than by the Fig. 33.-Cloud rms velocity versus cloud size for the same numerical cloud sample as in the previous figures, but with both rms velocity and cloud size measured in three-dimensional space, σv = σ v,3D / √ 3 and R cl . The thin solid line shows the average σv in logarithmic intervals of R cl , and the dashed line is the fit to the binned values, with a slope 0.37 ± 0.02 for cloud sizes R cl > 1.6 pc.
clouds self-gravity. This result is in contradiction to the claim that the velocity normalization of MCs scales with surface density ), based on clouds from the sample by Solomon et al. (1987), analyzed within rectangular maps of different sizes, rather than a fixed antenna-temperature threshold.
Besides being useful to derive the velocity normalization, the observed velocity-size relation may also provide a rather accurate estimate of the slope of the second order velocity structure function. In Figure 33 we show the relation between the velocity dispersion derived from the three-dimensional one, σ v,3D / √ 3, and the threedimensional cloud radius, R cl . The thick solid line is the least square fit to the average values of σ v in logarithmic bins of R cl shown by the thin solid line. The fits to the three-dimensional velocity-size relation for R cl > 1.6 pc (where the relation is well approximated by a power law) gives σ v = (1.19 ± 0.04) km s −1 (R cl /10 pc) 0.37±0.02 . (21) This relation is consistent with its two-dimensional counterpart (apart from the lower normalization due to the fact that R cl > R e on average), so the observable relation can be considered as a good estimate of the intrinsic three-dimensional one. Furthermore, the slope of the relation (21) is also consistent with the slope of the second order structure function, ζ(3)/2 ≈ 0.37, of the clouds from the simulation. Thus, we conclude that the observed velocity-size relation provides an estimate of the velocity scaling of MC turbulence, as long as it is based on MCs with reliable distance measurements and sufficient velocity resolution to detect the lower envelope of the relation.
The velocity-size relation (21) implies the following expression for the dynamical time of MCs as a function of the cloud three-dimensional radius, using the definition (13) of dynamical time adopted in §8: t dyn = 4.8 Myr (R cl /10 pc) 0.63 .
Using the result of §8 that t life /t dyn ≈ 4, our velocitysize relation implies that our largest MCs with sizes in

Mass-Size Relation
The mass-size relation is plotted in Figure 34, this time including all observed clouds above 100 M . The values of R e of the observed MCs have been divided by 10, to avoid overlap with the clouds from the simulation. Because of the imposed limit on the minimum cloud mass, the data is binned and fit only for R e > 2 pc for both the simulation and the observations. The resulting fits are M cl = (9.6 ± 0.3) × 10 3 M (R e /10 pc) 2.55±0.03 , (23) for the simulation, and M cl = (10.9 ± 0.7) × 10 3 M (R e /10 pc) 2.49±0.07 , (24) for the observations, so the slope of the mass-size relation from the simulation is consistent with the observations. The normalization of the mass-size relation depends on the threshold antenna temperature of the observational sample and the threshold density of the numerical sample. Of our MC catalogs described in §3, the ones with n H,min = 200 cm −3 have the mass-size normalization closest to that of the Outer-Galaxy MC sample by Heyer et al. (2001). This is the reason why all plots in this section are based on the highest-resolution catalog with that value of threshold density.
The total scatter in the relation is also similar between the observations and the simulation, if we account for the facts that the observational sample size is approximately twice as large as the numerical one, and that we have not added any observational uncertainty to the masses and sizes of the clouds from our simulation. Because of the similarity in both the slope and the total scatter, we can conclude that the mass-size relation resulting from SN-driven turbulence is consistent with that of real MCs from the Outer Galaxy Survey. A similar mass-size relation (though with a five times larger normalization) was also derived from the Galactic Ring Survey (Roman-Duval et al. 2010) and is implied by previous estimates of cloud fractal dimensions from various observational surveys (e.g. Elmegreen & Falgarone 1996;Sánchez et al. 2007) and from simulations of randomly driven turbulence (e.g. Kritsuk et al. 2007;Federrath et al. 2009).
Combining the mass-size relation (23) with the dynamical time expression in equation (22), and adopting the average value derived in §10.2 for the ratio between the two definitions of cloud radius, R e = 0.87R cl , we obtain the following expression for the average dynamical time of MCs as a function of cloud mass, hence the expression (18) for the average cloud lifetime as a function of cloud mass anticipated in §8.
10.5. Virial Parameter In §7, the virial parameter is computed with the threedimensional radius, R cl . Here, as in most MC studies including Heyer et al. (2009), we define an observable version of the virial parameter using the equivalent radius, R e , and thus refer to it as α vir,e . The dependence of α vir,e on M cl is fully determined by the velocity-size and mass-size relations presented above. Nevertheless, it is presented here as an alternative view of the comparison of the simulation with the observations, and to suggest an empirical calibration of the virial parameter of MCs.
The values of α vir,e are plotted versus M cl in Figure  35 for both the simulation (empty circles) and the observations (filled circles). The most striking feature of this plot is the very large scatter in the values of α vir,e at a fixed cloud mass, growing with decreasing cloud mass, which is a direct consequence of the large scatter in the velocity-size relation. As in that relation, the lower envelope for the observational data is limited by the smallest value of σ v to which the observations are sensitive, σ v ≈ 0.5 km s −1 . The dashed-dotted line in Figure 35 shows the minimum value of α vir,e as a function of M cl , computed by setting R e as a function of M cl using the average mass-size relation from the previous section and setting σ v = 0.5 km s −1 .
The virial parameter would be independent of mass if the velocity size relation were σ v ∼ R 1/2 e and the masssize relation were M cl ∼ R 2 e , the often assumed form of the Larson relations. However, the mass-size relation is such that the cloud surface density grows with mass, as shown in the previous section, causing the decrease of α vir,e with increasing M cl seen in Figure 35, only partly mitigated by the exponent of the velocity-size relation being a bit smaller than 0.5. The upper envelope in Figure 35 is even steeper than the average decrease of virial parameter with mass, as a consequence of the nearly flat upper envelope of the velocity-size relation. The dashed line in Figure 35 shows the expected upper envelope based on the mass-size relation and a maximum rms velocity of 2.5 km s −1 that is representative of some of the largest values in both the simulation and the observations.
Once we account for the minimum velocity dispersion in the Outer-Galaxy Survey, the relation of α vir,e with M cl and its scatter for the clouds from the simulation are consistent with those for the observed Outer-Galaxy MCs, as expected from the agreement found in the velocity-size and mass-size relations. This agreement suggests the possibility of an empirical calibration of the observed values of the virial parameter based on the results of our simulation. We have shown in §7 that the virial parameter computed with the radius R cl is α vir ≈ 1.2 (2E k /E g ). We have also shown in this section that, on the average, R cl ≈ R e /0.87, thus we derive this useful result for the relation between the observed virial parameter based on the equivalent cloud radius and the energy ratio: 2E k E g ≈ 0.96 α vir,e (assuming negligible saturation of the observed lines, such that the rms velocity can be assumed to be approximately mass weighted). Bertoldi & McKee (1992) modeled extensively the coefficient of the virial parameter of clumps, depending on the mass profile and the ellipticity of the clumps. Given their complex structure, MCs are not described by radial density profiles and ellipticity parameters as easily as compact clumps, so an empirical calibration based on simulations as proposed here is needed.

OVERVIEW OF RESULTS AND CONCLUSIONS
The main goal of this work is to test if ISM turbulence driven only by SN explosions can explain the turbulence observed within MCs. We have addressed this question with an AMR simulation representing an ISM volume of (250 pc) 3 and reaching a maximum resolution of 0.24 pc, with refinement based on density, density gradients and pressure gradients to resolve individual SN remnants and their interaction with the dense gas. We have studied the SN-driven turbulence over the whole volume and within individual clouds. We have compiled 12 different catalogs of MCs selected from the simulation using four different values of threshold density and three different spatial resolution. The properties of these clouds have been studied using tracer particles, hence taking advantage of the full resolution of the simulation. First we have presented cloud properties based on particle position, velocity and magnetic field values measured in three-dimensional space; then we have carried out a comparison with real MCs from the Outer Galaxy Survey by measuring projected quantities, such as the line-of-sight velocity dispersion and the equivalent radius. Our results are summarized in the following: 1. Near equipartition of total energies in the whole volume results in a distinct energy separation in the dense gas. While the overall ISM turbulence is trans-Alfvénic and mildly supersonic, the turbulence in the dense gas is highly supersonic and super-Alfvénic.
2. Approximately 11% of the total kinetic energy is transferred to the dense gas, even if most SN explosions occur at low densities. The dense gas has an average velocity dispersion of 8.5 km s −1 .
3. During the rapid expansion of a SN remnant, the velocity power spectrum may briefly show strong features at different scales. During most of the time the power spectrum is statistically relaxed and develops a power-law inertial range that scales with wavenumber as E(k) ∝ k −1.46 , between approximately 20 and 60 pc.
4. Unlike previous studies of compressible turbulence, the power spectrum of the compressive component of the velocity, E c (k) ∝ k −1.98 , is much steeper than that of the solenoidal component, E s (k) ∝ k −1.31 . The baroclinic effect is the best candidate to explain this result, as previous simulations of supersonic turbulence adopted equations of state where such effect was absent.
5. SN driving is not purely compressive. The curl of the forcing is the baroclinic term that is comparable to or larger than the divergence of the forcing (except in the unrealistic cases of a uniformdensity medium or a barotropic equation of state). As a result, the power in solenoidal modes exceeds that in compressive modes almost at any time and any wavenumbers (at the scale of 20 pc, E c /E s ∼ 0.2). Thus, isothermal simulations of turbulent fragmentation based on random solenoidal driving are a much better approximation of MC turbulence than isothermal simulations with purely compressive driving.
6. The time-averaged energy-injection scale of SNdriven turbulence is approximately 70 pc with our SN rate (may be larger with a smaller SN rate, or vice-versa), and oscillates in time between 50 and 100 pc.
7. The scaling exponents of the first and second order structure functions of velocity, relative to the third order one, in SN-driven turbulence are consistent with those found in supersonic turbulence driven by an idealized, random volume force, which supports the validity of turbulent fragmentation and star formation studies where the ISM turbulent cascade within MCs was modeled with such an idealized driving.
8. Based on a new scheme to compute densityweighted velocity structure functions, we obtain a third order exponent close to unity, as in the incompressible Kolmogorov's '4/5 law'. Thanks to the AMR method and to this density-weighting scheme, the structure functions probe the inertial range of MC turbulence down to a scale of 2-3 pc, while previous studies of SN-driven turbulence did not resolve an inertial range and only addressed the relative scaling.
9. The scaling of the velocity structure functions within individual MCs selected from the simulation is generally consistent with the scaling derived from MC observations. However, deviations are found for MCs directly affected by recent SN remnants.
10. The ratio of cloud virial parameter and kinetic to gravitational energy ratio is α vir /(2E k /E g ) = 1.2, independent of energy ratio and mass (see point 15 for the observable virial parameter). This structural property of MCs is not significantly affected by self-gravity during the duration of our simulation, where self-gravity is included in the final 11 Myr, corresponding to about two average cloud free-fall times. The ratio becomes very large only in a small fraction of clouds with α vir 10, where self-gravity causes the collapse of dense cores. Even in these star-forming clouds, there is no evidence of global cloud collapse. Self-gravity only causes a slight shift towards larger densities of the probability distribution of cloud densities.
11. The formation and dispersion times of MCs are of the order of two cloud dynamical times. Equivalently, the cloud lifetime, defined as the sum of formation and dispersion times, is approximately four cloud dynamical times. This is evidence indicating that SN-driven turbulence is responsible for cloud formation and dispersion, with little influence from self-gravity visible during the duration of our run. Future work, with longer runs, is needed to determine to what extent this remains true for longer periods of time.
12. The clouds have a mean magnetic field enhanced only by a factor of two relative to the mean magnetic field in the simulated volume The turbulence is super-Alfvénic for all clouds more massive than approximately 10 3 M .
13. The comparison with the MC sample from the Outer Galaxy Survey shows that clouds selected from the simulation have properties consistent with the observations, such as the mass and size distributions, the velocity-size and mass-size relations and the dependence of virial parameter on cloud mass. In our run, these properties, including the normalization of the velocity-size relation, are essentially the same for clouds selected either before or after the inclusion of gravity in the simulation; they are primarily the result of SN-driven turbulence, with only a minor contribution from selfgravity.
14. The normalization of the velocity-size relation does not depend on surface density in the simulation, nor in the observations. It is the same for MCs from the Outer Galaxy Survey as for those in the Galactic Ring Survey whose surface density is significantly higher.
15. The simulation provides a calibration of the observable virial parameter, α vir,e , based on the equivalent radius, which allows a derivation of the energy ratio from observational quantities, 2E k /E g ≈ 0.96 α vir,e .
Based on these results from the simulation and given its successful comparison with the Outer Galaxy Survey, we conclude that the SN-driven turbulence in our simulation is consistent with the observed MC turbulence during the duration of our experiment. Although other energy sources are present in galaxies, and local radiative and mechanical feedbacks also play a role in the dispersion of star-forming MCs, the origin, evolution and internal dynamics of MCs in our run are primarily the consequence of SN-driving, which is able to sustain turbulence at observed levels without help from those extra energy sources.
We are grateful to Christoph Federrath and the anonymous referee for providing useful comments and corrections. Computing resources for this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center, by PRACE through a Tier-0 award providing us access to the computing resource SuperMUC based in Germany at the Leibniz Supercomputing Center, and by the Port d'Informació Científica (PIC), Spain, maintained by a collaboration of the Institut de Física d'Altes Energies (IFAE) and the Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas (CIEMAT). PP acknowledges support by the ERC FP7-PEOPLE-2010-RG grant PIRG07-GA-2010-261359 and by the Spanish MINECO under project AYA2014-57134-P. TH is supported by a Sapere Aude Starting Grant from The Danish Council for Independent Research. Research at Centre for Star and Planet Formation was funded by the Danish National Research Foundation and the University of Copenhagen's programme of excellence.