# ISM properties in hydrodynamic galaxy simulations: Turbulence cascades, cloud formation, role of gravity and feedback

## Abstract

We study the properties of ISM substructure and turbulence in hydrodynamic (AMR) galaxy simulations with resolutions up to 0.8 pc and M. We analyse the power spectrum of the density distribution, and various components of the velocity field. We show that the disk thickness is about the average Jeans scale length, and is mainly regulated by gravitational instabilities. From this scale of energy injection, a turbulence cascade towards small-scale is observed, with almost isotropic small-scale motions. On scales larger than the disk thickness, density waves are observed, but there is also a full range of substructures with chaotic and strongly non-isotropic gas velocity dispersions. The power spectrum of vorticity in an LMC-sized model suggests that an inverse cascade of turbulence might be present, although energy input over a wide range of scales in the coupled gaseous+stellar fluid could also explain this quasi-2D regime on scales larger than the disk scale height. Similar regimes of gas turbulence are also found in massive high-redshift disks with high gas fractions. Disk properties and ISM turbulence appear to be mainly regulated by gravitational processes, both on large scales and inside dense clouds. Star formation feedback is however essential to maintain the ISM in a steady state by balancing a systematic gas dissipation into dense and small clumps. Our galaxy simulations employ a thermal model based on a barotropic Equation of State (EoS) aimed at modelling the equilibrium of gas between various heating and cooling processes. Denser gas is typically colder in this approach, which is shown to correctly reproduce the density structures of a star-forming, turbulent, unstable and cloudy ISM down to scales of a few parsecs.

###### keywords:

galaxies: evolution – galaxies: ISM – galaxies: star formation – galaxies: structure – ISM: kinematics and dynamics – ISM: structure## 1 Introduction

Observations of the interstellar medium (ISM) of nearby galaxies reveal a large variety of complex substructures, such as density waves, molecular clouds, filamentary structures along shocks, and supernova bubbles. ISM turbulence is a major driver of large-scale star formation: turbulent compression can increase the gas density and initiate star-forming instabilities (Elmegreen, 1993, 2002). This can happen in the spiral arms of modern disk galaxies (Klessen, 2001), and there is recent evidence that very active star formation relates to increased turbulence and associated instabilities in both primordial disk galaxies (Elmegreen et al., 2005; Förster Schreiber et al., 2006; Bournaud et al., 2007; Dekel et al., 2009) and starbursting mergers (Bournaud et al., 2008a; Teyssier, 2010). On the other hand, gas turbulence can also disrupt gaseous structures faster than they gravitationally collapse, which could potentially stabilize gas disks and quench the star formation activity of some galaxies (Martig et al., 2009; Dekel et al., 2009).

Nevertheless, the large-scale properties of ISM turbulence are still poorly understood and the ability of numerical simulations to realistically model ISM substructures remains largely unexplored. Hydrodynamic simulations of ISM turbulence and fragmentation in whole galaxy models have been performed (Wada et al., 2002; Wada & Norman, 2007; Tasker & Bryan, 2006; Tasker & Tan, 2009; Agertz et al., 2009a) but have not been compared to observational constraints yet, probably because of lacking the required level of realism, both in spatial resolution and in global dynamical modeling (no old stellar components and live dark halo were considered), leading to a restricted dynamical range of structure sizes. All these models are characterized by a probability distribution function (PDF) of the gas density that has a log-normal shape, resulting from the combination of various structures including density waves, dense gas clouds, supernovae bubbles. The levels of turbulent speed in these models seem realistic for both nearby and high-redshift galaxies (see Agertz et al., 2009a; Ceverino et al., 2010, respectively). Whether a realistic spectrum of structures of various scales is reproduced in modern hydrodynamic models remains however unknown.

Observations suggest that turbulent motions are initiated by gravitational instabilities on a relatively large scale (the Jeans scale, which is pc for wavenumber in nearby disk galaxies, Elmegreen et al. 2003) and that a turbulent cascade develops towards smaller scales. This turbulence should be three-dimensional since the Jeans scale length is expected to set the gas disk scale height. Turbulent motions could be triggered by gravity even inside molecular clouds (Field et al., 2009).

Observations of the power spectra of ISM emission from nearby galaxies are consistent with this. On scales smaller than pc, the power spectrum is steep, , as it is for 3D Kolomogorov turbulence, while on scales larger than pc, the power spectrum is flatter by about 1 in the slope, with a power-law exponent of . This has been observed in the Large Magellanic Cloud (LMC) (Elmegreen, Kim, & Staveley-Smith 2001; Block et al. 2010), and in NGC 1058 (Dutta et al. 2009). This tentatively suggests a two-dimentional regime for turbulent motions on large scales, because the break in the density distribution occurs at the expected disk scale height. The nature of the large-scale motions is uncertain, however, because the three-dimensional velocity field in individual galaxies is not observable. In the Small Magellanic Cloud, the power spectrum is steep everywhere (Stanimirovic et al. 2000), presumably because there is no thin disk on the line of sight.

In this paper, we explore the properties of ISM turbulence in simulations of isolated disk galaxies, performed with an adaptive grid hydrodynamic code (AMR) with a maximal resolution of 0.8 pc resolution. We use a “pseudo-cooling” approach to model the equilibrium of gas between various heating and cooling processes, based on a barotropic equation of state (EoS), and perform models with and without supernovae feedback. Our main model has a mass distribution and rotation curve representative for the LMC disk, so as to compare with observations from Block et al. (2010). However we do not try to reproduce the detailed individual structures of the real LMC, and we also analyze models with higher disk masses and gas fractions. We show that a large range of substructures spontaneously form in the ISM, with a power spectrum quite consistent with observations of the LMC from very large scales ( kpc) down to very small scales ( pc). A break in the power-law density power spectrum is observed at the gas disk scale height. Studying the velocity structure on various scales, we show that this corresponds to a transition from a 3D regime to a (quasi-)2D regime in the gas motion, as initially proposed from observations (Elmegreen et al. 2001). Even if large-scale forcings are probably important, the two-dimensional regime on large scales cannot be interpreted just as density waves and streaming motions in a rotating disk, and evidence for an inverse cascade of turbulence in a quasi-2D disk is proposed. The same properties are recovered in a model of a massive high-redshift disk galaxy, which has a high gas fraction (50%) and strong turbulence (), suggesting these properties are relatively universal for gaseous galactic disks.

Analyzing the vertical structure of the disk, compared to the properties of gas turbulence, we suggest that the turbulent motions at all scales are mainly powered by gravitational instabilities around the Jeans scale length, which sets the disk scale height. Feedback processes such as supernovae explosions do not significantly change the initial statistical properties of the ISM. They are nevertheless fundamental to maintain a realistic ISM structure over a large number of disk rotations. Supernovae bubbles do not appear to set the disk scale height, but they expand up to the disk scale height that is maintained by gravitationally-driven turbulence: this breaks apart the densest self-gravitating clouds to maintain a steady state cloud population.

## 2 Simulations

### 2.1 Simulation technique

We perform closed-box model of isolated galaxies, with physical parameters described in section 2.2. We use the AMR code RAMSES (Teyssier, 2002). Stars and dark matter are described with collisionless particles, evolved with a particle-mesh (PM) technique. The gaseous component is described on the same adaptive grid; hydrodynamic equations being solved with a second-order Godunov scheme.

Our simulations use a box size of 26 kpc. The coarse level of the AMR grid, , corresponds to a Cartesian grid, i.e. the minimal spatial resolution is pc where kpc is the cell size. The maximal refinement level in dense regions is set by , i.e. an effective resolution equivalent to and a maximal spatial resolution of pc. As we want to resolve turbulent motions in the whole disk, we use a refinement strategy that ensures that the mid-plane of the initial gas disk is entirely resolved with (pc), and with ( pc) at the initial average surface density. The coarsest resolution levels are used only in the outer stellar bulge and halo. The resolution increases towards a cell size of 0.8 pc in the inner disk and in dense substructures forming in the outer disk. Effective resolution maps are shown on Figure 3. The refinement is as follows: a cell of level is refined if its gas mass is larger than M and/or the number of particles (stars and dark matter) is larger than .

Only the gas density is computed in the smallest cells of the AMR grid with , the mass density of particles being restricted to in the PM scheme. This ensure a gravitational softening of at least 2 pc for particles. Treating the particle density at the level could result in a low number of particles per cell, and induce two-body relaxation, in regions that are refined because of a high gas density without necessarily having a high density of stars and/or dark matter. Only a few cells at each instant would actually contain less than 3 particles per cell, as measured in our simulation results, but we prefer avoiding any spurious effect.

Resolving the gas cooling and heating equation is numerically costly and would reduce the achievable resolution. In order to resolve small-scale structures and gas turbulence at the best possible resolution, we model the cooling and heating processes using an EoS that fixes the gas temperature as a function of its local atomic density . Here we use an EoS modeling the equilibrium of gas between the main heating (UV heating and stellar feedback) and cooling (atomic and dust cooling) processes (see also Teyssier, 2010). This EoS is a fit to the average equilibrium temperature of gas at one third solar metallicity: for densities cm, K. For densities cm, K. Densities cm correspond to gas in gravo-thermal equilibrium in the halo, at the Virial temperature K.

Gravitational instabilities are a major driver of turbulence in galactic disks, both by directly triggering gas motions and by fueling star formation that can also pump turbulent energy through feedback processes. This is known from observations (Elmegreen, 2002), and models similar to ours but scaled to massive spiral galaxies (Wada et al., 2002; Tasker & Tan, 2009; Agertz et al., 2009a). A fundamental aspect in our models is thus to avoid numerical instabilities and artificial fragmentation: a gas medium that should theoretically be stable, tends to be unstable in numerical models when the Jeans length is described by a small number of spatial resolution elements (Truelove et al., 1997). It is generally admitted that this effect is avoided if the Jeans length is resolved by at least a few elements (Machacek et al., 2003; Agertz et al., 2009a, e.g.). We then force the grid refinement to ensure that the thermal Jeans length of the gas is always resolved by at least 4 cells – and hence the actual thermal+turbulent Jeans length is resolved by an even larger number of cells). At the grid level, we impose a temperature threshold computed to keep the thermal Jeans length larger than . This temperature floor, added to the EoS at very high density ( M pc and above) can be considered as a sub-grid model for the unresolved turbulent motions at scales smaller than the maximal grid resolution. We can then consider that artificial fragmentation is generally avoided: simulations of galactic disk fragmentation imposing higher numbers of resolution elements per Jeans length do not show significant differences in the properties of gas stability and turbulence (e.g., Ceverino et al., 2010).

One of our models includes only disk self-gravity, but the other model also includes star formation and feedback. A threshold for star formation is implemented at a gas mass density atoms cm M pc. This choice ensures that stars form only in the densest gas structures formed in our models but that the density PDF remains correctly sampled around and above this threshold (see Figure 5). Above the threshold, the specific star formation rate (SFR per gas mass unit) is defined as the product of the efficiency and the local free-fall time. We use an efficiency of 10%, which was chosen to give a gas consumption timescale of 2 Gyr on the results of the first simulation without star formation. Note that lower resolution models usually need to combine a lower efficiency with a lower density threshold for star formation (see discussion in Teyssier et al. 2010).

We use the kinetic feedback model described by (Dubois & Teyssier, 2008). A 20% fraction of the mass of stars formed at a given timestep is assumed to be in stars that will explode into supernovae. 15% of the typical energy of each supernova, erg, is re-injected into the gas, assuming that the remaining energy is radiated away. The energy is re-injected in the form of radial velocity kicks around the supernova, in a gas bubble of radius 3 pc (see Dubois & Teyssier 2008 for details).

We use outflow boundary conditions for the fluid, and isolated boundary conditions for gravity

### 2.2 Model Parameters

Our simulations are not aimed at reproducing the exact morphology of the LMC, but at studying the properties of turbulence in a galactic disk with general properties similar to the LMC. To this aim, we initialize an exponential stellar disk containing particles, of total mass M radial scale-length of 1.5 kpc, and truncation radius of 3 kpc. A non-rotating bulge of M with radius 500 pc, made-up of particles, is added. The stellar velocity dispersions initially correspond to a stellar Toomre parameter of , so that stars are stable against axisymmetric instabilities but unstable to bar forming instabilities – a bar and a pair of spiral arms form spontaneously in our model.

The dark matter halo contains particles, has a total mass of M and follows a Burkert profile of core radius 3.5 kpc, truncated at 10 kpc. The initial gas disk is exponential, with a scale-length of 3 kpc and truncation radius of 5 kpc, an exponential scale-height of 50 pc truncated at 500 pc, and a gas mass of M. The gas disk is initially purely rotating, with no macroscopic velocity dispersion. The initial densities in the disk lie within the isothermal part of our EoS, so a thermal support of km s ( K) is present all over the initial disk, which can be considered as a model for the turbulent support of typically a few km s in real disks. This support ensures an initial Toomre parameter for the gas . Once dense structures form in the course of the simulation, the temperature decreases (following the EoS) and the thermal support is gradually replaced by a turbulent support. Outside the initial disk truncation, the AMR grid is initialized with a background density of atom per cm.

The rotation curve in this model is shown on Figure 2. It is broadly similar to the LMC rotation curve (Kim et al. 1998), with almost solid body rotation up to radii of at least 2 kpc, and gradually flattening with a circular velocity km s at large radii. We study the internal dynamics of the LMC; the interaction with the Milky Way and its hot gas halo is not modelled. Ram pressure was found to compress the leading edge of the HI and H disk, with little effect on the inner star-forming disk (Mastropietro et al. 2009; see also Tonnesen & Bryan 2009 on the effect of ram pressure on a cloudy ISM).

## 3 Results

The main simulation analyzed in Sections 3.1 and 3.2 is the LMC-sized model with stellar feedback. We later compare to the simulation without star formation and feedback in Section 3.3, to better highlight the role of feedback in the global structure of the ISM.

### 3.1 ISM structure and density power spectrum

We measured the velocity dispersion in the gas disk, and found that both the mass-weighted average value and the standard deviation rapidly increase in the first 200 Myr and do not significantly evolve after Myr, indicating that a steady state in the turbulent regime has been reached. We thus mostly analyze our models around =250 Myr.

Spectral analysis of the gas distribution was performed on a sub-box of size 13 kpc (the outer regions of the full simulation box contain only the outer regions of the dark matter halo). Throughout the paper, the wavenumber is for a wavelength of 13 kpc, which results in a wavenumber unit pc.

Figure 3 shows the gas density distribution for this simulation at Myr. The response to a stellar bar, and some long spiral arms are seen, together with more chaotic structures such as shorter arms, shocks, dense clouds, bubbles, etc (Fig. 4). We fit a sech vertical profile to the gas at radii smaller than 5 kpc, integrated over the plane, and find a scale height of 207 pc. The scale height is relatively constant, with moderate flaring, and varies from 187 pc inside the central kpc to 226 pc for kpc.

The density PDF, integrated over the central disk with a radius of 5 kpc and height of 2 kpc, is shown on Figure 5. This PDF has a log-normal shape, truncated at quite high density ( cm) which is enabled by our high spatial and mass resolution. If we write the probability distribution for the density

then the Gaussian dispersion is .

The power spectrum of the face-on gas column density is shown on Figure 6 at the same instant and later ones. No significant time evolution is found. The power spectrum shows a double power-law with a break at wavelengths of about 150 pc, which is about the measured gas disk scale-height. The power law for large structures has a slope of and the slope for small structures is . The slope steepens for the smallest scales below 5-10 pc, which is discussed later in section 3.3.

The simulation without star formation and feedback shows a globally similar density power spectrum (see Sect. 3.3 and Fig. 14) with a double power-law shape with slopes of on small scales and on large scales. The transition occurring around the gas disk scale-height. Only the smallest scales below 10 pc significantly differ from the model with star formation and feedback (see section 3.3).

Our simulations produce a density structure that has about the power spectrum observed for the ISM in the LMC and other nearby galaxies. The power spectrum shows a break between a power law and a power law, with the transition occurring at scales around the gas disk scale-height, as usually speculated in observations (e.g., Elmegreen et al., 2001; Dutta et al., 2009). In the following, we study the velocity structure to understand the nature of the motions in the two regimes of the density distribution.

### 3.2 Velocity structure

total | ||

Velocities (km s) |

total | ||

Velocities (km s) |

#### Velocity fields and power spectra

Maps of the in-plane radial velocity component and perpendicular component are shown on Figure 8 for the same snapshot as Figures 3 and 4, all at Myr. Another instant ( Myr) is shown on Figure 9 for the gas density and Figure 10 for the velocity maps. The velocity maps are shown for the total velocity components, and for the low wavenumbers ( pc) and the high wavenumbers ( pc) separately. This separates the two regimes identified on the density power spectrum, where the transition occurs at or pc. These maps were built by zeroing the low- or high- components in the result of a Fourier transform of the velocity field, and recovering the velocities through an inverse Fourier transform.

The power spectrum of the three velocity components , and is shown on Figure 7. A single power law is found for the in-plane components. The perpendicular velocity follows the same power law on small scales, but its power spectrum flattens on large scales, with a transition at about the same scale length (150 pc) as in the density structure spectrum.

These differences in velocity power spectra correspond to higher power for the long-wavelength in-plane components than the long-wavelength perpendicular components of velocity. There are large in-plane motions at long wavelengths from spiral and bar-driven gas flows and relatively little perpendicular motions at long wavelengths to accompany them.

#### A 3D cascade of turbulence on small scales

On scales smaller than the gas disk scale height (100-200 pc), the ISM density has the density power spectrum expected for a 3D cascade of turbulence. The three-dimensional nature of the associated motions is confirmed by the high-wavenumber velocity maps, which show about the same amplitude for the in-plane motions and the perpendicular one. Typical average values of on various scales are given in Table 1. The small-scale motions are globally isotropic^{1}^{2}

These results are overall consistent with a 3D cascade regime on the small-scale part of the density spectrum, initiated at scales around 150 pc, which is the gas disk scale height and the average Jeans length. This is naturally the case if the most gravitationally unstable scale, namely the Jeans length, corresponds to the disk thickness. Indeed, the gas disk scale height is expected to be set by the gaseous Jeans length in nearby disk galaxies (Elmegreen et al., 2001; Dutta et al., 2009); this is also the case in high-redshift disk galaxies, which are more turbulent, have larger Jeans lengths, and have thicker gas disks at the same time (Elmegreen & Elmegreen, 2006; Genzel et al., 2006; Förster Schreiber et al., 2006; Bournaud et al., 2008b). We checked this in our simulations by computing a mass-weighted average of the velocity dispersion (8 km s), the average surface density , and then an “effective” Jeans length pc (locally, the Jeans length can significantly differ from this average value, as the density and velocity dispersion vary). This is consistent with the disk scale height being set by the Jeans length.

Time | 254 Myr | 268 Myr |
---|---|---|

small scales () | 1.13 | 1.07 |

large scales () | 0.16 | 0.14 |

total velocities (all ) | 0.92 | 0.89 |

#### A 2D inverse cascade on large scales?

The low-wavenumber (large-scale) structures have a flatter density power spectrum. They also have different kinematic properties. The large-scale component of the velocity fields is dominated by in-plane motions. Perpendicular motions are much weaker. The low-wavenumber component of shows only modest peaks in Figure 8, and these peaks are generally found in low-density regions (they generally correspond to gaseous fountains above the disk plan, rather than vertical motions in the mid-plane).

The motions corresponding to the large-scale regime of the density power spectrum are quasi-2D, highly anisotropic (Table 1). This is not because the disk rotation dominates these motions: this applies to as much as , and the ratio is almost unchanged when we remove the lowest wavenumbers tracing global rotation. The power spectra of the three velocity components (Fig. 7) actually shows that the perpendicular velocity is much lower than the in-plane components for all wavelengths larger than the disk scale-height.

The slope of the large-scale branch of the density power spectrum and the quasi-2D nature of the associated motions suggest that we are observing something like a 2D inverse cascade of turbulence (Kraichnan, 1967). We nevertheless have to explore the nature of these motions more accurately, to ensure that these properties are not a conspiracy from density waves and associated gas flows, which would then not be turbulent motions.

The large-scale components of the in-plane motions (Figs 8 and 10) show a central mode, which is the characteristic response to a bar+arms system of density waves (e.g., Combes & Gerin, 1985; Athanassoula, 2000). This response is however quite asymmetric, and does not dominate the large-scale motions outside of the central kpc. Overall, outside the central kpc, chaotic motions dominate the large-scale velocity components (of course superimposed on the disk rotation pattern). In particular the comparison of the gas density and velocity maps shows large vortices around dense gas clumps: for instance, the densest clump seen to the left of Figure 3 is clearly associated with an in-plane eddy, with positive and negative peaks of surrounding the position of this dense clump. Such signatures are reminiscent of a Rosby Wave instability (Lovelace et al., 1999; Varnière & Tagger, 2006), which is observed in pure 2D disks but also in thick disks or 3D systems (Meheut et al., 2010) and could thus take place in our quasi-2D system. Such motions would be quasi-2D turbulence, which is more than a large-scale flow in a density wave, but which may be induced by a density wave of other large-scale forcings.

The vorticity map of the large-scale velocity component (), after subtraction of the rotation pattern, is shown on Figure 11. It shows a number of vortices with various sizes and strengths throughout the disk. The power spectrum of the enstrophy is shown in Figure 12. Enstrophy is the integral of the square of the vorticity: . The enstrophy has a power-law power spectrum throughout the entire range of 2D scales, like the in-plane motions shown before. The slope of this power spectrum is -1, just as expected for 2D turbulence, as found in two-dimensional experiments (e.g., Petersen et al., 2006, 2007). This enstrophy spectrum means that the quasi-2D part of the density power spectrum is made-up of vortices with size and strength distributed as expected for 2D-turbulence: this is unlikely to simply result from non-turbulent flows along density waves.

The properties of the gaseous motions on scales larger than the disk scale height thus contain a component from turbulence in a quasi-2D medium. The origin of this turbulence could be a combination of long-range forcing from spiral waves and other global disturbances, and an inverse cascade from smaller scales, where energy is put into the ISM by gravitational instabilities on the Jeans length and stellar feedback. Energy injection at the Jeans length is consistent with our observation of vorticity around the gravitating clumps and with the quasi-2D power spectrum on scales larger than the Jeans length. Wada et al. (2002) already mentioned that the properties of turbulent motions in their galaxy models were consistent with an inverse cascade of turbulence. Their models were purely two-dimensional, while we here suggest a similar inverse cascade in a three-dimensional gas disk, which is thin but has a well resolved scale height.

The large-scale motions probably cannot simply result from a single inverse cascade with energy input only at the Jeans scale. Indeed, the growth rate of gravitational perturbations in a disk decreases only as where is the scale of the perturbations, so that long-range forcing at scales larger than the typical Jeans scale length remains significant. At least, the gravitational coupling of stellar and gaseous components and the bi-component stability (Jog & Solomon, 1984) should result in a large range of unstable scales, ranging typically from the gaseous Jeans scale to the stellar Jeans scale (Elmegreen, 1995). There should then be a range of unstable scales injecting energy into the inverse turbulence cascade. While the expected properties of 2D turbulence, such as the enstrophy spectrum, are recovered in our simulations, one can note that the slope of the density power spectrum in this regime is slightly higher than expected from pure two-dimensional turbulence (a power law, instead of the exponent expected for purely two-dimensional and uncompressible turbulence). The non-zero thickness of the gas disk and the relative importance of long-range forcing in a gaseous+stellar disk might cause this.

The role of long-range gravitational forcing and the coupling of the gaseous and stellar components also appear by comparing our simulations with the AMR models in Tasker & Tan (2009). The stellar disk is modeled with a rigid analytical potential in their simulations. Gas clumps form everywhere in their disk with a relatively peaked size distribution and a constant separation of about one Jeans scale. The Jeans length is also a characteristic scale in our model: it sets the disk scale height and gas clouds along spiral arms are also often separated by (Fig. 4), but gaseous structures overall form over a wider range of sizes in our models. The inclusion of both stellar and gaseous gravity on large scales appears crucial to realistically reproduce the full spatial range of the main ISM structures from long spiral arms and large vortices to smaller clumps and shocks, even if ISM structures on small scales can be realistically reproduced without stellar gravity (e.g., de Avillez & Breitschwerdt, 2007).

We have analyzed above the LMC-sized simulation with star formation and stellar feedback. In this simulation, the properties of ISM turbulence seem mainly driven by gaseous and stellar gravity. However, the comparison to the model without stellar feedback below will highlight a fundamental role of feedback in maintaining a steady state in the density distribution of the ISM.

### 3.3 Role of star formation feedback in gas disk properties

The “gravity-only” simulation (without star formation and feedback) shows a relatively similar distribution of gaseous structures, and a relatively similar density power spectrum compared to the model with feedback (see Figures 13 and 14), at least at Myr and at scales larger than 5-10 pc. Besides suggesting that our earlier results do not crucially depend on a specific feedback scheme, it also indicates that the ISM structuring and the pumping of turbulent energy into the turbulent cascades results mainly from gravitational processes and is not primarily driven by supernovae explosions. The disk scale-height and the disk substructures are unchanged in this gravity-only model.

Feedback processes nevertheless play a fundamental role in maintaining this gravity-driven turbulence distribution. We can indeed note already at Myr in Figure 14 that the ISM power spectrum shows an excess of very small-scale structures (below 5-10 pc) and the gas density map shows the associated small and dense gas clumps. Over longer timescales the power spectrum becomes less realistic as the large scales are gradually emptied and the bump at very small scale increases: gas dissipates its energy and accumulates in tiny dense bullets from where it is never removed (Fig. 13).

The role of feedback appears when comparing the density power spectra in our two models: feedback disrupts the dense smallest-scale structures, and supernovae bubbles expand up to the disk scale height. Density maps show supernovae-driven bubbles, all of which are smaller than or comparable to the disk scale-height. But this does not mean that the size of supernovae bubbles regulates the scale height: the measured scale-height and the break in the density power spectrum are not changed by feedback; the scale height is mostly regulated by gravitational processes, and is about the Jeans length. But feedback disrupts structures on the smallest scales and expands them up to the Jeans length, where gravitational processes can initiate a new cycle of 3D turbulent cascade towards small scales and 2D inverse cascade towards large scales.

Overall, the gas disk structure and its scale height are largely regulated by gravitational instabilities, but gas structures and energy dissipation on the smallest scales are also regulated by star formation feedback. The break in the face-on density power spectrum and the scale height measured on edge-on projections both correspond to the average Jeans length. 2D and 3D motions are initiated by instabilities around this scale length and large-scale forcings. Feedback processes prevent the otherwise inevitable accumulation of gas in tiny and dense bullets, disrupt small-scale structures and refill the turbulent cascades initiated at the Jeans length. This cycle maintains the ISM density distribution in a steady state.

Observations have already suggested that turbulent motions are primarily driven by gravitational instabilities such as spiral arms and molecular clumps (Elmegreen et al., 2003). Based on the analysis of the density distribution, these authors pointed out the role of star formation feedback in breaking apart small dense structures to preserve a steady state. Agertz et al. (2009a) also proposed from simulations that the main driver of ISM turbulence could be gravitational forcing, with nevertheless a significant contribution from star formation (see also Dib, Bell & Burkert 2006).

Gravity-only models without star formation and feedback initially produce a realistic density distribution in the ISM over a couple of disk rotation periods. But on longer timescales and over large numbers of disk rotations, the inclusion of feedback in numerical models appears necessary to maintain a realistic density distribution. However, this applies only to high-resolution simulations where gas cooling is modeled down to K. Lower-resolution SPH simulations, and more generally models where the ISM is only modeled as a smooth stable gas at K or more (e.g., the stabilizing EoS by Springel & Hernquist 2003 and Robertson et al. 2004), would not be affected by catastrophic gas dissipation if they do not include feedback sources, as they do not resolve the turbulent and unstable nature of the star-forming ISM.

### 3.4 Extension to massive and gas-rich galaxies

The results described so far were for a low-mass disk galaxy, scaled to the LMC properties. To explore whether the highlighted properties can extend to more massive disk galaxies and/or disk galaxies with different gas fractions, we have applied the same analysis to a high-redshift galaxy model described in Delaye et al. (2010 in prep.). This is a massive galaxy (baryonic half-mass radius of 10 kpc and circular velocity of 280 km s) with a 50% gas fraction, modeled with a 2 pc resolution. As typical star-forming galaxies at , it has a high gas fraction, a very unstable disk with high turbulent speeds (a few tens of km s) and giant complexes of star formation in a thick ( kpc) gas disk (see observations by Elmegreen & Elmegreen (2006); Elmegreen et al. (2007); Genzel et al. (2006) and simulations by Semelin & Combes (2005); Bournaud et al. (2007); Agertz et al. (2009b); Ceverino et al. (2010)).

The density power spectrum is shown on Figure 15. It is consistent with a double power law, with the break occurring at wavelength around 1 kpc which is the disk scale height and the typical Jeans length in such galaxies. We measured a total ratio of velocity dispersions of , with on large scales and on small scales. The large-scale motions are still significantly non-isotropic, even if the low wavenumber velocities are somewhat more isotropic than in our LMC-sized model, which could be caused by the disk being thicker. The properties in the density distribution and chaotic motions found in our LMC-sized models thus seem to apply, at least qualitatively, to galaxies with different masses, rotation curves and gas fractions.

### 3.5 Resolution requirements

Our fiducial LMC-sized model with star formation and feedback has been re-done with , i.e. a spatial resolution of 13 pc, which means that the disk scale height is resolved with 10 resolution elements, instead of 100-200 resolution elements at the maximal refinement level . The density threshold for star formation was reduced by a factor of 80 to keep the resulting star formation rate similar.

The global density distribution of this reduced-resolution model is not very different from the initial case, except that the highest density structures in the PDF are removed; otherwise the PDF is only moderately shifted towards lower densities and has about the same dispersion (Fig. 5). However, the velocity structure is significantly different. The initial model had almost isotropic velocity dispersion (the smallest scales begin the most isotropic ones). Interestingly, the reduced-resolution model has much less isotropic motions, with (at Myr and over all wavenumbers; the ratio was 0.9 in the high-resolution run). We thus note that a very high resolution is required to resolve the isotropy or non-isotropy of turbulent motions in disk simulations, and that a few (10) resolution elements per disk scale height are not sufficient. Our interpretation is as follows: vertical motions can be initiated or amplified when instabilities occur off the midplane at heights 0 or 0, when the Jeans length becomes locally smaller than the disk scale height (which happens in dense regions or low-velocity dispersion regions). This requires to resolve, locally, several Jeans length per scale height. Then, each Jeans length itself needs to be resolved by a least a few resolution elements: 4 resolution elements per thermal Jeans length is a strict minimum imposed in our model and our EoS generally results in 5-10 elements per thermal Jeans length, and the number of resolution per thermal+turbulent Jeans length would be even higher. Thus, capturing instabilities at and the associated triggering of vertical motions can require a few tens of resolution elements across the gas disk scale height. When the disk thickness is marginally resolved, the non-isotropy of gas motions can be artificially increased at small scales. In this case, each gas clump is mostly spinning around its minor axis (see Agertz et al. 2009a) generally aligned with the spin axis of the entire galaxy (but see Tasker & Tan 2009).

## 4 Summary and conclusions

### 4.1 Properties of ISM turbulence, role of gravity and feedback

This paper has presented an LMC-sized models for comparison with observations of the LMC (Block et al. 2010) and a high-redshift clumpy disk model. The LMC-sized model does not match all individual structures observed in the LMC – this is not the initial purpose and the LMC has many low-density ionised filaments, and other substructures in low-density regions, that are out of reach of our mass resolution.

The model is successful in reproducing the double power law shape of the ISM density power spectrum observed in the LMC (Elmegreen et al. 2001 and Block et al. 2010) in a system with a similar mass, gas fraction and rotation curve. We analyzed the 3D velocity structure of the models, which could not be done in observations, to understand the origin of this density substructures.

Our results indicate that a dual turbulence cascade could take place in our simulations, and by extension in the real ISM. The scale height of the gas disk is set by the Jeans scale length, which is the most gravitationally unstable scale. Gravitational instabilities at/around this scale length inject energy in turbulent motions. These motions follow a classical 3D cascade of isotropic turbulence towards smaller scales, on which they form a hierarchy of gas clouds, GMCs and substructures therein. Turbulent motions on scales down to and inside GMCs thus appear to be from self-gravity. Nevertheless, stellar feedback plays a major role in disrupting dense structures on the smallest scales and in maintaining the turbulence cascade in a steady state.

On scales larger than the gas disk scale height, very anisotropic turbulent motions with a shallower density power spectrum are observed. Long-range gravitational forces are probably not negligible, but the vorticity and enstrophy properties on these large scales are as expected for (quasi-)2D turbulence. The large scale motions are not limited to global disk rotation and streaming flows along density waves. This inverse cascade appears to be initiated mainly by the same gravitational instabilities at the Jeans length as the direct 3D cascade, because it is present with and without stellar feedback.

While models without self-gravity do not find a characteristic scale for energy injection into ISM turbulence (Joung & Mac Low, 2006), we identify a major role of gravitational processes around the Jeans length. Large-scale forcing and stellar feedback at small scales are important processes in regulating the ISM properties too, but the Jeans length is the scale at which chaotic gas motions change from quasi-2D to 3D isotropic turbulence, and at which a break occurs in the density power spectrum. Log-normal gas density PDFs have been found in models of supernovae-driven ISM turbulence without self-gravity (e.g. de Avillez & Mac Low, 2002), but without having the observed power spectrum reproduced. The power spectrum analysis in our models and in LMC observations rather suggests that the energy injection in the ISM is largely from gravitational processes (such as gravitational instabilities and inward mass accretion, Elmegreen & Burkert, 2010; Klessen & Hennebelle, 2010) with a regulating role of feedback on small scales.

We have shown that similar properties for ISM substructures and turbulence, at least qualitatively, are found in models of galaxies with higher masses and higher gas fractions. This includes the extreme case of high-redshift () galaxies with thick, highly turbulent and very clumpy gaseous disks. An interesting implication for observers is that the isotropy of the gas velocity dispersion is strongly dependent on the scale, and hence on the resolution of an observation.

### 4.2 Modeling ISM substructures in galaxy simulations

We have explored ISM modeling with a barotropic equation of state (EoS) corresponding to gas in thermal equilibrium between a standard radiation field with the main atomic and fine-structure cooling processes. We have shown that this EoS employed in an AMR hydrodynamic code with a quasi-Lagrangian refinement strategy reproduces a realistic power spectrum for the ISM density distribution. It does not form a truly multiphase medium, since the pressure is only a function of density. It nevertheless a clumpy ISM, with cold and dense clouds embedded in a warmer diffuse phase, and the resulting density distribution is consistent with that of the real, multiphase ISM. This technique is efficient to reach high spatial and mass resolutions: the computation of cooling and heating rates is not numerically costly, but the absence of out-of-equilibrium, colder or warmer gas in this model prevents the timescale to become extremely short, which would typically happen in parsec-scale models with complete cooling calculations.

Hydrodynamic models with cooling calculations produce a log-normal PDF for the ISM (Wada & Norman, 2007; Tasker & Tan, 2009; Agertz et al., 2009a), which is retrieved in our model over a larger dynamical range. The PDF alone does not guarantee that the correct mass fraction is distributed in structures of various sizes and densities. Our model with a mass and rotation curve similar to the LMC matches the density spectrum observed by Block et al. (2010) in the LMC. While all thermal processes are not explicitely computed, this at least means that the ISM modeled this way has realistic density distribution, with a correct number of gas structures of various sizes and a correct mass fraction therein, from the global galaxy-sized scales down to 5-10 pc, i.e. molecular clouds and their main substructures.

The ISM substructures are mostly powered by gravity-driven turbulence and instabilities around the Jeans scale length. Star formation feedback is however an essential ingredient, not only by structuring the ISM into shells locally, but also by providing an extra source of energy is needed to disrupt the smallest and densest structures that tend to be steadily filled by gas dissipation into compact clumps. Gravity, hydrodynamics, and a realistic thermal model are initially enough to produce a realistic structure distribution in the ISM: feedback seems relatively negligible over a couple of disk rotations, but becomes important over longer timescales to balance systematic gas dissipation. The source of energy that disrupt the structures at very small scales and refills the 3D cascade and the inverse 2D cascade at the Jeans scale is, in our model is supernova feedback. Other sources such as HII regions, stellar winds, and radiation pressure from young massive stars (Murray et al., 2010; Krumholz & Dekel, 2010) could also contribute to maintain the ISM density distribution in a steady state. The global density structure of the ISM remains relatively unchanged in simulations where the disk scale height is marginally resolved, but modelling accurately the quasi-isotropic turbulent motions on small scales is found to require a very large number (at least a few tens) of resolution elements per disk scale height.

## Acknowledgments

This work was granted access to the HPC resources of CCRT and CINES under the allocation 2010-GEN2192 made by GENCI (Grand Equipement National de Calcul Intensif). We are grateful to John Scalo and Sebastien Fromang for suggestions, and to Oscar Agertz, Daniel Ceverino, Françoise Combes and Elizabeth Tasker for discussions and comments on the manuscript. DLB is indebted to Mr F. Titi as well as the AVENG Group of Companies for their sponsorship. DLB thanks Roger Jardine, Kim Heller and the AVENG Board of Trustees. IP acknowledges support from the Mexican foundation Conacyt.

### Footnotes

- The slightly higher values for the vertical velocity component correspond to low-density gas flows ejected above the disk plane by supernovae winds; the simulation without feedback has values of even closer to for high wavenumbers.
- Note that is slightly larger than , so the in-plane motions are not perfectly isotropic. This is expected for instance if the gas tends to follow epicyclic motions, for which the ratio of the radial to tangential velocity dispersions is .

### References

- Agertz, O., Lake, G., Teyssier, R., Moore, B., Mayer, L., & Romeo, A.B. 2009a, MNRAS, 392, 294
- Agertz, O., Teyssier, R., & Moore, B. 2009b MNRAS, 397, L64
- Athanassoula, E. 2000, in Stars, Gas and Dust in Galaxies: Exploring the Links, ASP Conference Series, Vol. 221 p. 243, Danielle Alloin, Knut Olsen & Gaspar Galaz Eds.
- de Avillez, M. A., & Mac Low, M.-M. 2002, ApJ, 581, 1047
- de Avillez, M. A., & Breitschwerdt, D. 2007, ApJL, 665, L35
- Block D. L., Puerari I., Elmegreen B. G. & Bournaud F., ApJL, 718, L1
- Bournaud, F., Elmegreen, B. G., & Elmegreen, D. M. 2007, ApJ, 670, 237
- Bournaud, F., Duc, P.-A., & Emsellem, E. 2008, MNRAS, 389, L8
- Bournaud, F., et al. 2008b, A&A, 486, 741
- Delaye, L., et al. in preparation
- Ceverino, D., Dekel, A., & Bournaud, F. 2010, MNRAS, 440
- Combes, F., & Gerin, M. 1985, A&A, 150, 327
- Dubois, Y., & Teyssier, R. 2008, A&A, 477, 79
- Dutta, P., Begum, A., Bharadwaj, S., & Chengalur, J. N. 2009, MNRAS, 397, L60
- Dekel, A., Sari, R., & Ceverino, D. 2009, ApJ, 703, 785
- Dib, S., Bell, E., & Burkert, A. 2006, ApJ, 638, 797
- Elmegreen, B. G. 1993, ApJl, 419, L29
- Elmegreen, B. G. 1995, MNRAS, 275, 944
- Elmegreen, B. G., Kim, S., & Staveley-Smith, L. 2001, ApJ, 548, 749
- Elmegreen, B. G. 2002, ApJ, 577, 206
- Elmegreen, B. G., Elmegreen, D. M., & Leitner, S. N. 2003, ApJ, 590, 271
- Elmegreen, B. G., Elmegreen, D. M., Vollbach, D. R., Foster, E. R., & Ferguson, T. E. 2005, ApJ, 634, 101
- Elmegreen, B. G., & Elmegreen, D. M. 2006, ApJ, 650, 644
- Elmegreen, D. M., Elmegreen, B. G., Ravindranath, S., & Coe, D. A. 2007, ApJ, 658, 763
- Elmegreen, B. G., & Burkert, A. 2010, ApJ, 712, 294
- Ferland, G. J., Korista, K. T., Verner, D. A., Ferguson, J. W., Kingdon, J. B., & Verner, E. M. 1998, PASP, 110, 761
- Field, G., Blackman, E., Keto, E. 2009, MNRAS submitted, arXiv:0904.4077
- Förster Schreiber, N. M., et al. 2006, ApJ, 645, 1062
- Genzel, R., et al. 2006, Nature, 442, 786
- Jog, C. J., & Solomon, P. M. 1984, ApJ, 276, 114
- Joung, M. K. R., & Mac Low, M.-M. 2006, ApJ, 653, 1266
- Kim, S., Staveley-Smith, L., Dopita, M. A., Freeman, K. C., Sault, R. J., Kesteven, M. J., McConnell, D. 1998, ApJ, 503, 674
- Klessen, R. S. 2001, ApJ, 556, 837
- Klessen, R. S., & Hennebelle, P. 2010, A&A submitted (arXiv:0912.0288)
- Kraichnan, R.H. 1967, Phys. Fluids, 10, 1417
- Krumholz & Dekel 2010 MNRAS in press, arXiv:1001.0765
- Lovelace, R. V. E., Li, H., Colgate, S. A., & Nelson, A. F. 1999, ApJ, 513, 805
- Machacek, M. E., Bryan, G. L., & Abel, T. 2003, MNRAS, 338, 273
- Martig, M., Bournaud, F., Teyssier, R., & Dekel, A. 2009, ApJ, 707, 250
- Mastropietro, C., Burkert, A., & Moore, B. 2009, MNRAS, 399, 2004
- Meheut, H., Casse, F., Varniere, P., & Tagger, M. 2010, A&A, 516, A31
- Murray, N., Quataert, E., & Thompson, T. A. 2010, ApJ, 709, 191
- Petersen, M.R., K. Julien, & J.B. Weiss 2006, Physics of Fluids, 18, 26601
- Petersen, M.R., K. Julien, & G.R. Stewart: 2007, ApJ, 658, 1236
- Robertson, B., Yoshida, N., Springel, V., & Hernquist, L. 2004, ApJ, 606, 32
- Semelin, B., & Combes, F. 2005, A&A, 441, 55
- Springel, V., & Hernquist, L. 2003, MNRAS, 339, 289
- Stanimirovic, S., Staveley-Smith, L., van der Hulst, J. M., Bontekoe, T.J. R., Kester, D. J. M., & Jones, P. A. 2000, MNRAS, 315, 791
- Tasker, E.J., & Bryan, G.L. 2006, ApJ, 641, 878
- Tasker, E.J., & Tan, J.C. 2009, ApJ, 700, 358
- Teyssier, R. 2002, A&A, 385, 337
- Teyssier, R., Chapon, D., & Bournaud, F. 2010, ApJL submitted, arxiv:1006.4757
- Tonnesen, S., & Bryan, G. L. 2009, ApJ, 694, 789
- Truelove, J. K., Klein, R. I., McKee, C. F., Holliman, J. H., II, Howell, L. H., & Greenough, J. A. 1997, ApJl, 489, L179
- Varnière, P., & Tagger, M. 2006, A&A, 446, L13
- Wada, K., Meurer, G., & Norman, C.A. 2002, ApJ, 577, 197
- Wada, K., & Norman, C. A. 2007, ApJ, 660, 276