MASS TRANSPORT AND TURBULENCE IN GRAVITATIONALLY UNSTABLE DISK GALAXIES. I. THE CASE OF PURE SELF-GRAVITY

Until the past few years, most theoretical work on galaxy evolution has focused on reproducing observed correlations in the bulk properties of galaxies, such as the stellar mass-halo mass relation (e.g., Behroozi et al. 2013; Moster et al. 2013), the stellar mass–star formation relation (e.g., Daddi et al. 2007; Elbaz et al. 2007; Noeske et al. 2007), and the mass–metallicity relation (e.g., Pilyugin et al. 2004; Tremonti et al. 2004). In such correlations, galaxies are treated as single points, and theoretical models have been largely content to treat them as such. In these models, the behavior of galaxies is almost entirely dictated by a balance between cosmological accretion and feedback (e.g., Bouché et al. 2010; Lilly et al. 2013; Forbes et al. 2014b; Mitra et al. 2015), and for this reason most theoretical attention has been focused on stellar feedback as the dominant driver of galaxy evolution.

The recent availability of large samples of spatially resolved maps of gas, metals, and star formation in nearby galaxies has opened up a new frontier in the study of galaxies: explaining the radial distribution of these quantities. Observational studies conducted to date have turned up puzzling facts that demand explanation, which cannot obviously be explained simply by a local balance of accretion versus feedback at all locations in a galaxy. Instead, these observations point to a potentially significant role for radial redistribution of material as a key process for understanding the evolution and growth of galaxies over cosmological timescales.

The first observational puzzle concerns the radial distribution of gas and star formation in present-day disk galaxies. Observations show that the neutral gas in such galaxies is distributed with an approximately universal exponential profile, with typical gas scale lengths of $\sim 0.5\;{r}_{25}$, where r₂₅ is the optical radius of the galaxy (Regan et al. 2001; Schruba et al. 2011; Bigiel & Blitz 2012). The gas is dominated by ${{\rm{H}}}_{2}$ inside $\sim 0.4\;{r}_{25}$, while H i predominates at larger radii (e.g., Leroy et al. 2008). The radial distribution of star formation is, to first approximation, simply linearly proportional to the ${{\rm{H}}}_{2}$ distribution (Bigiel et al. 2008; Schruba et al. 2011; Leroy et al. 2013), so star formation is much more radially concentrated than the total neutral gas, and predominantly occurs inside $\sim 0.5\;{r}_{25}.$ In such regions, the time required to convert all the available gas to stars is $\sim 2\;{\rm{Gyr}}$, much less than a Hubble time (Bigiel et al. 2008; Leroy et al. 2013). Thus the inner parts of galaxies will become gas-depleted and cease star formation unless fresh gas is supplied at an equal rate, which for L_* galaxies is a few ${M}_{\odot }$ yr⁻¹ (Kennicutt & Evans 2012). While some disk galaxies, including the Milky Way, do have inner gas holes and quenched central star formation, such systems appear to be the exception rather than the rule among the local disk galaxy population (Bigiel & Blitz 2012).

How the fresh gas responsible for the lack of quenched centers can be explained is an unsolved problem. While condensation of gas from the hot halo (Marinacci et al. 2010; Fraternali et al. 2013; Hobbs et al. 2013) will supply gas in a centrally concentrated manner, direct accretion of cold gas from the intergalactic medium (e.g., Kereš et al. 2005; Dekel et al. 2009a) primarily deposits gas at large galactocentric radii, far from the actively star-forming regions that need to be refueled. This can lead to build-up of gas at large radii (c.f. Dutton 2012). Gas recycling from evolved stellar populations can potentially provide some of the required mass (Leitner & Kravtsov 2011) at small galactocentric radii, but, in addition to uncertainties about whether this channel provides enough mass and is consistent with various chemical evolution constraints, the radial distribution of recycled gas has not been explored. In summary, the lack of quenched, gas-depleted centers in local disk galaxies seems to require some form of gas redistribution within the disk, for which gravitational instability is an obvious candidate.

A second surprising observation regarding the radial distribution of gas concerns gas velocity dispersions. H i lines in disk galaxies show super-thermal velocity dispersions of $\sim 10$ km s⁻¹ (e.g., van Zee & Bryant 1999; Petric & Rupen 2007; Tamburro et al. 2009; Ianjamasimanana et al. 2012, 2015). Velocity dispersions are highest toward the centers of galaxies, but they decline only shallowly with radius, and remain superthermal even well outside r₂₅. While star formation feedback alone appears able to drive the observed velocity dispersion in inner disks, the same is not true outside r₂₅, where the star formation rate drops precipitously but the velocity dispersion does not (Tamburro et al. 2009; Ianjamasimanana et al. 2015). Dwarf galaxies represent an extreme in this regard: due to their very low rates of star formation, supernova feedback cannot plausibly provide enough energy to drive the observed turbulence anywhere within them (Stilp et al. 2013). Nor do magneto-rotational or thermal instability appear to be sufficient (Kim et al. 2003; Piontek & Ostriker 2004, 2005, 2007; Yang et al. 2007). Again, gravitational instability is an obvious candidate to drive the turbulence in the weakly star-forming portions of galaxies (although see Elmegreen & Hunter 2015, who show that gravitational instability is not globally important in very low-mass dwarfs).

A final surprising observation concerns the metallicity distributions of galaxies. Observations of disk galaxies show that they have slightly negative gas-phase metallicity gradients inside r₂₅ (e.g., Vila-Costas & Edmunds 1992; Pilyugin et al. 2004; Henry et al. 2010; Balser et al. 2011; Ho et al. 2015), and nearly flat distributions of metallicity outside r₂₅ (Bresolin et al. 2009, 2012; Werk et al. 2011). Dwarf galaxies show no metallicity gradients at all (Croxall et al. 2009). While the gradients in inner disks might be explicable simply via a balance between inflow, outflow, and star formation (e.g., see Portinari & Chiosi 2000 and Spitoni & Matteucci 2011 versus Ho et al. 2015), this is not the case for the flat gradients in outer disks. In these strongly gas-dominated regions, the metallicity should simply scale with the stellar mass fraction, regardless of the presence of either inflow or outflow. While the stellar mass fraction changes sharply with galactocentric radius outside r₂₅, the metallicity does not, a finding that seems extremely difficult to explain without invoking some form of metal redistribution. As with the previous two results, gravitational instability seems a candidate redistribution mechanism worthy of further exploration.

Given the likely importance of gravitational instability as a driver of galaxy evolution, it is not surprising that there have been a number of theoretical studies devoted to it. Some of the earliest were one-dimensional (1D) models by Lin & Pringle (1987a, 1987b), Olivier et al. (1991), and Ferguson & Clarke (2001), who argued that the observed exponential distribution of stars in galactic disks could only be understood if gas undergoes significant viscous transport on a timescale comparable to the timescale over which the stellar disk forms. These models assumed a fixed dimensionless viscosity, but Krumholz & Burkert (2010; generalizing earlier work by Bertin & Lodato 1999 in the context of Keplerian disks with no stars) showed that the viscosity and thus the rates of mass and angular momentum transport could be computed self-consistently by balancing energy loss against turbulence generation by gravitational instability.

Building on this work, (Forbes et al. 2012, 2014a) argue that the observed structure of star forming disks is a natural result of an equilibrium between gravitational instability, accretion, and star formation. Cosmological accretion brings gas to the disk outskirts, where the gas is in general stable to collapse. Over the course of a few rotation periods, gravitational instability drives torques that tend to move gas to smaller galactocentric radii. The inward flow of gas feeds active ongoing star formation in the inner disk, such that the total star formation in the inner disk is ultimately modulated by the accretion of gas at large galactocentric radii. In these models the bulk of the disk is at all times near ${Q}_{{\rm{total}}}\sim 1.$

Our motivation for three-dimensional (3D) simulations is that, while the 1D models are instructive and useful for exploring parameter space quickly, they have substantial limitations. The gas and stars both contribute to the self-gravity of the disk, so the degree of gravitational instability of both components must be considered (Lin & Shu 1966; Jog & Solomon 1984a, 1984b; Bertin & Romeo 1988; Wang & Silk 1994; Romeo & Wiegert 2011), in the form of the Toomre (1964) stability criterion. While this is possible in 1D models, it requires strong assumptions about the vertical and azimuthal structure of the disk, and about the response of each component. Moreover, material at a given radius is only allowed to move inward or outward at one bulk velocity. Properly modeling the full nonlinear multi-component gravitational instability process requires 3D modeling of both the gas and stars embedded in a realistic galactic gravitational potential field.

Recently, a new class of high-resolution numerical galaxy formation simulations have arisen that begin to probe the density and resolution scales necessary to resolve star forming clouds. These models include a variety of simplifying assumptions, including an analytic stellar background potential, both with and without imposed spiral arms (Dobbs et al. 2006; Tasker & Tan 2009; Tasker 2011; Smith et al. 2014; Dobbs et al. 2015; Tasker et al. 2015). Others resolve the gas to very high spatial resolution at the cost of capturing only a fraction of a galactic dynamical time, or of simulating a galaxy much smaller than the Milky Way (Bournaud et al. 2010; Renaud et al. 2013). Still others take a similar approach to ours (Agertz et al. 2009; Hopkins et al. 2012; Fujimoto et al. 2014a, 2014b; Behrendt et al. 2015), including runs based on the same AGORA (Kim et al. 2014) initial conditions that we use (Agertz et al. 2015).

In this paper we make use of simulations of the gas, stars, and dark matter in idealized isolated models of star forming disk galaxies. We focus on the spatial and mass scale similar to the Milky Way and study three different disks with varying gas fractions. We expand on previous work in a number of ways:

• 1.  
  Using our high numerical resolution—${\rm{\Delta }}x\sim 20\;\mathrm{pc}$ at the maximum refinement level, corresponding to a threshold density for star formation, ${n}_{{\rm{thresh}}}\sim 50\;{\mathrm{cm}}^{-3}$—we are able to comfortably resolve the global properties of star forming regions in our disk models (although we do not resolve the internal properties of star forming clouds), allowing us to resolve an ISM with multiple distinct phases.
• 2.  
  Our simulations are run for several galactic dynamical times. This long time baseline allows the instability to reach a saturated, fully chaotic state. As we will show, the onset of violent gravitational instability in our disk models is accompanied by substantial transient ringing behavior characteristic of the smooth unstructured initial conditions, and we run for long enough for this initial state to be forgotten.
• 3.  
  We perform a detailed modeling of the Toomre Q parameter, taking into account contributions of the gas and stars, including both the thermal and kinetic gas pressure, and accounting for the nonzero thickness of the gas and stellar disks.
• 4.  
  We measure the radial and temporal variation in the rate of radial mass transport at high spatial and temporal resolution, allowing us to track the detailed dynamics of the flow of matter in the disk and measure a mean flow of matter through the disk.
• 5.  
  The initial conditions data, full simulation code, simulation outputs, and analysis code are publicly available.³

Our simulations are quite similar to those of Agertz et al. (2015). Their simulations are based on nearly identical initial conditions but were run with a factor of two higher spatial resolution than we use. On the other hand, we simulate galaxies over much longer baselines (600 Myr versus 140 Myr), which is essential to allow the gravitational instability to settle reach full nonlinear saturation and thereby to allow an accurate measurement of the rates of the typical mass transport rate due to gravitational instability. Agertz et al. (2015) do not consider the issue of transport in their study.

In addition, our simulations can be compared to the numerical experiments of Behrendt et al. (2015). They focused on the transition from the linear to nonlinear phase of gravitational instability in the context of an analytic treatment of the Toomre instability. We see similar morphology in our simulations—an initially smooth disk that breaks up into axis-symmetric rings which then further segment and collapse into bound clumps. Their simulations were run with an isothermal equation of state and did not include a stellar disk and furthermore did not consider radial mass transport.

This paper is part of a series. In the companion paper (N. J. Goldbaum et al. 2015, in preparation) we perform simulations including stellar feedback, and compare these to the pure-gravity simulations presented here. A previously published paper in this series (Petit et al. 2015) makes use of one of the suite of simulations we describe here in order to study metal transport by gravitational instability, and thus we will not examine that topic further in this work.

The plan for the remainder of this paper is as follows. In Section 2 we describe the code, subgrid physics models, and initial conditions. In Section 3, we describe in general terms the time evolution of our simulations, particularly with respect to the dynamical evolution of the gravitational instability process. In Section 4 we quantitatively measure the effect of the gravitational instability via the time evolution of the Q parameter, as well as the radial mass flux, measure the mean mass flux rate through our simulated disks, and discuss the velocity structure of the gaseous disk. Finally in Section 5 we review our findings and discuss implications for future work.

The initial conditions were generated for the AGORA project (Kim et al. 2014) using the makegalaxy code (Springel et al. 2005). Briefly, makegalaxy transforms an input halo mass (M₂₀₀), stellar disk mass (M_*), gas fraction (f_g), halo spin parameter (λ), and halo concentration parameter (f_c) into particle initial conditions. These input parameters uniquely set the rotation curve, ${V}_{C}(R)$, and radial scale length, h_R, via the analytic theory of Mo et al. (1998).

The dark matter and stars are represented in the galaxy model using collisionless particles and are initialized by stochastically drawing from analytic distribution functions. The dark matter positions are initialized to follow a Hernquist (1990) distribution, which closely matches the more commonly used Navarro et al. (1996) fitting formula, but is more analytically tractable. The star particle positions are initialized following

$$\begin{eqnarray}&&{\rho }_{*}(R,z)=\displaystyle \frac{{M}_{*}}{4\pi {h}_{z}{h}_{R}^{2}}{{\rm{sech}}}^{2}\left(\displaystyle \frac{z}{2{h}_{z}}\right)\mathrm{exp}\left(-\displaystyle \frac{R}{{h}_{R}}\right),\end{eqnarray} \tag{ 1 }$$

where M_* is the mass of the stellar disk, h_z is the vertical scale height, z is the vertical coordinate, and $R=\sqrt{{x}^{2}+{y}^{2}}$ is the cylindrical radial coordinate. We consider gas fractions ${f}_{g}=0.1$, 0.2, and 0.4, which we refer to as the low gas fraction (LGF), fiducial, and high gas fraction (HGF) runs, respectively (Table 1). Once the particle positions are computed, the velocities are populated using a distribution function that depends only on the local orbital energy E and vertical component of the angular momentum L_z. The velocity distribution function is assumed to be axissymetric, making it strraightforward to locally solve the Jeans equation making use of the known density distribution from above (see Springel et al. (2005) for more details).

In addition to the stellar disk component, we also include a bulge of mass M_B. The bulge is initialized following a Hernquist profile:

$$\begin{eqnarray}&&{\rho }_{{\rm{B}}}(r)=\displaystyle \frac{{M}_{{\rm{B}}}}{2\pi }\displaystyle \frac{b}{r{(r+{h}_{{\rm{B}}})}^{3}},\end{eqnarray} \tag{ 2 }$$

where h_B is a free parameter. The bulge-to-total mass ratio, ${M}_{{\rm{B}}}/({M}_{{\rm{B}}}+{M}_{*})$, is 0.1 in all three runs.

We include ${N}_{{\rm{H}}}={10}^{7}$ halo particles, ${N}_{{\rm{D}}}={10}^{7}$ particles in the stellar disk, and ${N}_{{\rm{B}}}=1.25\times {10}^{6}$ particles in the stellar bulge. All particles in each population have uniform masses. For the halo population, the mass of each particle is ${m}_{{\rm{H}}}=1.3\times {10}^{5}\;{M}_{\odot }.$ For the bulge and disk population the particle masses are ${m}_{{\rm{D}}}={m}_{{\rm{B}}}=3.4\times {10}^{3}\;{M}_{\odot }.$

Since makegalaxy produces initial conditions formatted for the gadget smoothed particle hydroynamics (SPH) code, special care must be taken to initialize the gas onto Enzo's AMR grid structure. Rather than interpolating from the initial conditions for the SPH particles, we instead initialize the gas density following an analytic density profile:

$$\begin{eqnarray}&&\rho (R,z)=\displaystyle \frac{{M}_{{\rm{G}}}}{4\pi {h}_{R}^{2}{h}_{z}}\mathrm{exp}\left(\displaystyle \frac{-R}{{h}_{R}}\right)\mathrm{exp}\left(\displaystyle \frac{-| z| }{{h}_{z}}\right),\end{eqnarray} \tag{ 3 }$$

where ${M}_{{\rm{G}}}={f}_{g}{M}_{*}$ is the mass of the gaseous disk. The gas velocities are initially axisymmetric, with the radial profile set according to the circular velocity curve written out by makegalaxy.

This departs somewhat from the procedure used by makegalaxy to generate SPH initial conditions, where the gas density distribution is assumed to be exponential in the radial direction but the vertical extent is determined via an iterative relaxation process to ensure the gas disk is initially in equilibrium. This relaxation process assumes that the gas pressure is moderated by a subgrid effective equation of state model, which we do not employ in this study. In practice, the gas in our simulations initially experiences a phase of violent collapse and relaxes into a quasi-equilbrium state after the gravitational instability has fully developed. This means the precise state of the gas in the initial conditions is not terribly important—the simulation "forgets" the gas initial conditions. Furthermore, this choice also allows us to make an apples-to-apples comparison with future isolated galaxy simulations that make use of the public AGORA isolated galaxy initial conditions.

Rather than making use of an analytic dark matter and stellar potential (see, e.g., Dobbs et al. 2006; Tasker & Tan 2009), we employ a live dark matter halo and stellar disk. This allows us to follow the active response of the stars and dark matter to the collapse of the gaseous disk. The particle dynamics are implemented in Enzo using a standard particle-mesh scheme (Hockney & Eastwood 1988). Particle positions and velocities are updated according to the local gravitational acceleration using a drift-kick-drift scheme.

Hydrodynamics are captured using the piecewise parabolic method (PPM; Colella & Woodward 1984). PPM hydrodynamics allows us to accurately capture strong shocks in a few computational zones while also maintaining second-order spatial accuracy. While the PPM method is formally second-order accurate in space and time, Enzo's adaptive timestepping scheme is only first-order accurate, so these calculations are only first order accurate in time. In addition, after experiencing intermittent instability using second-order interpolation, we opted to degrade to first-order interpolation at AMR level boundaries. This means that we are formally first-order accurate in space, although these inaccuracies should only show up at level boundaries.

Since the ∼200 km s⁻¹ circular velocity of the galaxy necessitates strongly supersonic flows in the galactic disk, we make use of the dual energy formalism implemented in the Enzo code (Bryan et al. 1995). In standard PPM hydrodynamics, the internal energy is not tracked—instead it is derived from the total energy and velocity. When the kinetic energy is much larger than the internal energy, this can lead to spurious temperature fluctuations due to floating point round-off error. To avoid this, Enzo tracks a separate internal energy field that in effect provides extended precision to the total energy field. This allows us to safely resolve the thermal physics of the gas in a global simulation of a disk galaxy without worrying about spurious temperature fluctutions due to advection errors.

Our simulations make full use of the AMR capabilities of the Enzo code. Our refinement strategy is to focus computational effort on the highest refinement level, while still adequately resolving the dark matter halo so our estimate of the gravitational potential in the neighborhood of the galactic disk is accurate. Since the dark matter halo is spatially extended compared to the disk, with the furthest dark matter particles living 500 kpc away from the center of the galaxy, this necessitates a hierarchically nested AMR structure.

To fully encompass the dark matter halo, we employ a cubic simulation box with a width of 1.3 Mpc, and resolve the root grid with 64³ cells. In addition to the static root grid, we impose 5 additional levels of statically refined regions, enclosing volumes that are successively smaller by a factor of 8.

In addition to the static refinement, we allow for an additional five levels of adaptive refinement. To keep the particles properly resolved at all times, we refine a cell if the total mass in particles within the cell exceeds $1.7\times {10}^{6}$ ${M}_{\odot }$, or approximately 10 halo particles. This choice produces nested grid hierarchies in the region of the simulation dominated by the dark matter halo where the grid hierarchy is determined adaptively.

To keep the gaseous disk resolved at all times, we refine a cell if the total mass in gas within the cell exceeds $2.2\times {10}^{4}$ ${M}_{\odot }$. Finally, to avoid artificial fragmentation, we ensure the Jeans length, ${\lambda }_{{\rm{J}}}=\sqrt{(\pi {c}_{s}^{2})/(G\rho )}$, is locally resolved at all times by at least 32 cells, comfortably satisfying the Truelove et al. (1998) criterion. To keep the Jeans length resolved after collapse has reached the maximum refinement level, we employ a pressure floor such that the Jeans length is resolved by at least four cells on the maximum refinement level. To avoid contaminating the hydrodynamics in situations where we can use more resolution, the pressure floor is only employed for cells on the maximum refinement level. The pressure floor is implemented as a source of "extra" pressure above the thermal pressure—an improvement over the Jeans stabilization routines in the current stable version of Enzo where the pressure floor is implemented by increasing the temperature. The new pressure floor implementation will be included in the next stable release of Enzo.

The Enzo AMR hierarchy is fully adaptive in space and time, so regions that no longer satisfy the refinement criteria described above will no longer be refined. We do not include an explicit de-refinement criteria, instead allowing the refinement criteria described above to fully control the AMR hierarchy.

The particle refinement criteria we use are not as stringent as the gas refinement criteria. If the gas were not present in the simulation, the star refinement criteria alone would cause the stellar disk to be refine to AMR level 7 (160 pc) resolution in the outskirts of the disk and AMR level 8 (80 pc resolution in the inner disk and bulge. Since the combined contribution of the gas and stars determine the final resolution of the stars, in general the stellar disk is resolved to $20-40\;\;\mathrm{pc}$ resolution. In particular, the resolution is higher near dense clumps and filamentary structure in the gas. Early in the simulation, before the gas has had a chance to be consumed by star formation, the resolution is also generally higher.

Near the beginning of all three simulations, there are $3\times {10}^{7}$ computational zones in total, with about $1.5\times {10}^{7}$ zones on the maximum refinement level. Toward the end of the fiducial and high gas fraction simulations, after much of the gas disk has been consumed the number of zones on the maximum refinement level decreases to $\sim 5\times {10}^{6}$ zones. We attribute this decline to the collapse of gas into dense clumps and the consumption of gas by star formation, leaving much of the volume occupied by relatively rarefied gas.

Since the gaseous disks in our simulations are unstable to fragmentation and collapse, it is necessary to include a subgrid model that converts dense, collapsing gas into newly formed star particles. Rather than using the built-in "standard" star formation prescription in the Enzo code (Cen & Ostriker 1992; The Enzo Collaboration et al. 2014), which is tuned for lower-resolution cosmological simulations, we make use of a new star formation prescription. To ease comparison with future papers using the AGORA initial conditions, the prescription is based on the suggested star formation model for the AGORA project (Kim et al. 2014).

Briefly, the model assumes that the star formation rate density in any cell is a function only of the gas density in that cell, according to the following formula:

$$\begin{eqnarray}\displaystyle \frac{d{\rho }_{*}}{{dt}}=\left\{\begin{array}{lr}{f}_{*}\displaystyle \frac{\rho }{{t}_{{\rm{ff}}}} & :\rho \gt {\rho }_{{\rm{thresh}}}\\ 0 & :\rho \leqslant {\rho }_{{\rm{thresh}}},\end{array}\right.\end{eqnarray} \tag{ 4 }$$

where f_* is the star formation efficiency, ${t}_{{\rm{ff}}}=\sqrt{3\pi /32G\rho }$ is the local dynamical time, and ${\rho }_{{\rm{thresh}}}=\mu \;{m}_{{\rm{h}}}\;{n}_{{\rm{thresh}}}$ is the threshold density for star formation. Here μ is the fixed mean molecular weight, m_h is the mass of a hydrogen atom, and ${n}_{{\rm{thresh}}}=50\;{\mathrm{cm}}^{-3}$ is the threshold number density. Note that since we do not advect species fractions for hydrogen and helium ionization states, $\mu =1.4$ is constant over the full simulation box. The value we choose is typical of the bulk of the atomic ISM. In all of our simulations, we use ${f}_{*}=1\%,$ in accordance with the observed low star formation efficiency universally observed in star forming regions at a range of density and size scales (Krumholz & Tan 2007; Krumholz et al. 2012).

The star formation threshold is chosen such that the Jeans length for gas on the maximum refinement level remains resolved by at least four cells until the refinement reaches the maximum level, assuming a temperature floor of 100 K, approximately the minimum temperature we see in our simulations. Once gas reaches the maximum allowed refinement level, the gas is no longer well resolved, and we allow the gas to convert itself into stars and apply a pressure floor to prevent artificial fragmentation. Since we only barely resolve the typical densities for the formation of molecular gas, we do not employ a subgrid model to track the molecular gas fraction, instead assuming that gas above the threshold density is fully molecular. This follows Guedes et al. (2011), who found that the precise value chosen for ${\rho }_{{\rm{thresh}}}$ was not important for determining the star forming properties of a disk galaxy in a cosmological zoom-in simulation, so long as the threshold density is high enough that star formation occurs primarily in dense clumps of gas rather than in the bulk of the ISM.

Practically speaking, Equation (4) cannot be solved by spawning star particles at each timestep for all cells with densities above the threshold density. Given a typical timestep of 1000 years, and assuming a cell containing gas at the threshold density, we would expect to be spawning star particles with masses of only 0.01 ${M}_{\odot }$. Rather than spawning star particles in every cell that exceeds the density threshold, we instead impose a minimum star particle mass, and form stars stochastically. In our scheme, a cell will form a star particle of mass ${m}_{{\rm{sf}}}=300\;{M}_{\odot }$ with probability

$$\begin{eqnarray}{P}_{*}=\left\{\begin{array}{lr}{f}_{*}\displaystyle \frac{\rho {\rm{\Delta }}{x}^{3}}{{m}_{{\rm{sf}}}}\displaystyle \frac{{dt}}{{t}_{{\rm{ff}}}} & :\rho \gt {\rho }_{{\rm{thresh}}}\\ 0 & :\rho \leqslant {\rho }_{{\rm{thresh}}},\end{array}\right.\end{eqnarray} \tag{ 5 }$$

where dt is the simulation timestep on the maximum refinement level. By stochastically creating according to Equation (5), the time-averaged value of the star formation rate in each cell is still given by Equation (4), while creating a reasonable number of particles.

In Figure 1, we compare the star formation rate inferred from the ages of star particles present at the end of each simulation with the expected star formation rate measured by evaluating Equation (4) at 25 Myr intervals. We see that there is good agreement between the measured and predicted star formation rate, with a maximum deviation of ∼7%. The implementation of this star formation algorithm has been made publicly available and will included in the next stable release of Enzo.

Zoom In Zoom Out Reset image size

To model the thermal physics of the gas in our simulations, we make use of the Grackle cooling library⁴ (Kim et al. 2014; The Enzo Collaboration et al. 2014). The Grackle cooling and heating routines were adapted from the chemical and thermal physics implementation included in the Enzo code. Relative to Enzo's implementation of thermal physics, Grackle adds novel capabilities that we make use of in this study. In addition, this choice eases comparison with future simulations performed in other codes as well as simulations performed for the AGORA comparison.

Grackle includes a primordial cooling routine based on tabulated cooling rates as a function of density and temperature, avoiding the need to evolve separate fields for each hydrogen and helium ionization state. In addition to cooling by primordial species, we also include cooling due to metal line emission. These cooling rates are inferred from tabulated rates output by the CLOUDY code.

The disk is initialized to include a metal color field that follows the initial gas distribution. All gas zones in the disk are initialized with a solar metal fraction. The metal field is passively advected along with the gas. We include a dynamical metal color rather than assuming a fixed metal fraction in anticipation of modeling the production of metals in simulations with supernova feedback.

To properly model the thermal physics in a Milky Way-like ISM, we include a prescription for heating due to electrons released from dust grains by the photoelectric effect. This is implemented as a constant heating rate of $8.5\times {10}^{-26}\ {\rm{erg}}\ {{\rm{s}}}^{-1}$ per Hydrogen atom applied uniformly throughout the simulation box for gas below ${10}^{4.3}$ K (Tasker 2011). This rate is chosen to match the expected heating rate assuming a UV background consistent with the solar neighborhood value (Draine 2011).

In order to analyze our simulations, we construct an extensive post-processing pipeline using a combination of the yt toolkit (Turk et al. 2011) and a set of custom python analysis scripts that we have made publicly available on Bitbucket.⁵ Details of the analysis pipeline, and how we reconstruct various quantities of interest from the raw simulation outputs, are described in Appendix A. In particular, Appendix A gives our formal definitions of quantities that are non-trivial to calculate, including velocity dispersions (Appendix A.4), effective sound speeds (Appendix A.5), and Toomre Q parameters (Appendix A.6).

Since our simulated galaxies do not include a prescription for star formation feedback, the primary driver for the dynamical evolution of our model galaxies is gravitational instability. While the initial conditions for our simulations are formally stable, cooling allows the gas to quickly lose hydrostatic support, leading to catastrophic collapse in the vertical direction. In addition, our shearing self-gravitating disks are susceptible to the Toomre (1964) instability.

This story can be inferred by inspecting Figure 4, where we plot the time evolution of the azimuthal average of ${Q}_{{\rm{total}}}.$ The initially stable disk quickly becomes unstable (e.g., the regions that show up in red in the bottom half of each subplot). The instability leads to a radially expanding wave of collapsing gas. Soon after, the gas collects in gravitationally bound clumps, which proceed to migrate through the disk, both radially inward and outward. The regions inside the gravitationally bound clouds are formally unstable according to a local Toomre analysis due to their very high surface densities. The interclump regions reach relatively high values of ${Q}_{{\rm{total}}}$ so these regions are formally stable to collapse.

Zoom In Zoom Out Reset image size

This process plays out in all three simulations, albeit with varying collapse speeds and degrees of violence. The high gas fraction case initially develops filaments, but by ∼200 Myr of evolution these have broken up into giant clumps everywhere in the disk. As we lower the gas fraction, the transition to the clump-dominated phase takes longer, and the clumps themselves become smaller. Nonetheless, all three simulations reach similar clump-dominated regimes, which appears to represent the fully saturated state of the gravitational instability.

We can see what is driving the evolution of ${Q}_{{\rm{total}}}$ by inspecting Figure 5. We show the evolution in ${Q}_{{\rm{total}}}$ (top row) along with the quantities that determine ${Q}_{{\rm{total}}}:$ the gas and stellar surface density and effective sound speed, as well as the epicyclic frequency. We see that both the epicyclic frequency and the stellar surface density and velocity dispersion show little variation over the course of the simulations.

Zoom In Zoom Out Reset image size

On the other hand, the gas surface density and effective sound speed show significant variation. The effective sound speed tends to increase in time and the gas surface density tends to decrease. Both of these effects drive a secular increase in ${Q}_{{\rm{gas}}}$ over the course of the simulation. In turn, this increase in ${Q}_{{\rm{gas}}}$ leads to the secular increase in ${Q}_{{\rm{total}}}$ Eventually, once the gas supply is exhausted, ${Q}_{{\rm{gas}}}\gg 1$, so ${Q}_{{\rm{total}}}$ approaches ${Q}_{*}.$

Here we examine the radial flow of material in the disks of our simulated galaxy models. We measure if there is any net flow of material either from the galactic center outward or from the outskirts inward. In addition, we examine the detailed radial and time dependence in the gas mass flux.

In Figure 6, we present our measurements of the gas mass flux as a function of radius and time. All three simulations exhibit similar overall behavior. Initially, the models exhibit significant ringing as the disks initially collapse and settle down. In this stage the mass flux is dominated by rings of material experiencing alternating bands of inward and outward flow. After the initial collapse phase, the gas collects in gravitationally bound clouds. For the rest of the simulation, the mass flux rate is primarily determined by the inward and outward flow of the gravitationally bound clumps. As spiral arms develop, the gravitational potential in the disk begins to develop non-axisymmetric components that tend to drive the gas clumps both radially inward and outward.

Zoom In Zoom Out Reset image size

Along with the detailed variation in the mass flux, we would also like to know if there is any net mass flux once the disk has settled down into a quasi-equilibrium state. To answer this question, we make use of the time averaging algorithm described in Appendix A.9. In this way we only consider the mass fluxes measured above the blue lines in Figure 6, which delineate the approximate time at which the disk has settled into a steady state.

This averaging procedure results in the time-averaged mass-flux measurements in the top panel of Figure 7. For each simulation, we plot the time-averaged mass flux as a function of radius. We find that there is a net inward flux of gas at most radii for all three model galaxies. The flux scales very roughly with the gas fraction, with a typical mass flux of −0.3, −1.0, and −2.0 ${M}_{\odot }$ yr⁻¹ mass flux in the LGF, fiducial, and HGF cases, respectively.

Zoom In Zoom Out Reset image size

Although the flow is inward at most radii in the disk, we note that the outer edges of the disks instead show an outward mass flux. (In Figure 7, this feature is outside the plotted region for the HGF case, but is visible for the Fiducial and LGF cases.) Comparing with Figure 4, we see that this outward mass flux appears roughly at the boundary between the gravitationally unstable region and the still-stable outer region of the disk. Such a reversal of the mass flow direction at the transition between the stable and unstable regions of a gravitationally unstable disk is consistent with the predictions of Forbes et al. (2014a). However, we caution that these outer regions have been averaged over for the shortest times, and so the exact value of the typical outward flux is uncertain.

We might expect that the inward flow of gas is sufficient to supply the star formation in the inner regions of our model galaxies. To see if this is the case, we plot the time-average of the radially accumulated star formation rate in the bottom panel of Figure 7. The mass flux in the top panel is sufficient to supply the star formation in the bottom panel only if the mass flux is greater than the cumulative star formation rate at any given radius. Due to the very high star formation rates in these simulations with no feedback (c.f. Figure 1), we find that the inward flow of gas is insufficient to fuel the star formation rate in these simulations.

However, we note that this is a result of the unphysically high star formation rates that our simulations exhibit due to the lack of star formation feedback. If we instead consider the star formation rate that would be expected given the gas surface density distribution combined with observed star formation rates (e.g., Leroy et al. 2013), we reach the opposite conclusion: our inflow rates are sufficient to fuel star formation at observed levels. Indeed, our fiducial, Milky Way-like simulation produces mass transport at a rate of $\sim 1$ ${M}_{\odot }$ yr⁻¹, which is roughly the observed star formation rate in the Milky Way (e.g., Chomiuk & Povich 2011). The question of whether inflows and star formation can be matched simultaneously in a simulation including feedback we defer to the companion paper.

Here we examine the detailed gas velocity structure in our simulated galaxies. We are particularly interested in the radial dependence in the effective sound speed, the anisotropy in the velocity dispersion, and the relative contribution of turbulent motions and thermal sound speed to the effective sound speed. In Figure 8 we present the time evolution of the azimuthally averaged effective sound speed in each of our simulations. In all three simulations, gravitational instability alone is able to drive substantial turbulent velocity dispersions in the inner disk. The outer, gravitationally stable disk is characterized by much lower velocity dispersions. We also note that the typical velocity dispersion scales with the initial gas fraction. This may be interesting for studies of high-redshift gas disks, which are known to exhibit substantial turbulent velocity dispersions.

Zoom In Zoom Out Reset image size

To obtain a more complete picture of the velocity structure in our simulations, we look in more detail at the time-averaged velocity structure in the fiducial simulation. We perform the same sort of time averages described in Appendix A.9, focusing on the velocity dispersion, sound speed, effective sound speed, and velocity anisotropy.

We present the results of this time averaging in Figure 9. In the top panel, we plot the average effective sound speed as a function of radius, along with the 1σ dispersion in the effective sound speed at any given radius. The effective speed is typically of order ∼20 km s⁻¹, with a gradual radial decline comparable to what is observed in H i maps of nearby galactic disks (see the discussion and references in Section 1.1). This is perhaps surprising given that feedback processes are often invoked as an energy source for turbulence. Instead, we find generically that gravitational instability alone is sufficient to drive substantial turbulent velocity dispersions, more than enough to match observed gas velocity dispersions in resolved observations of galactic gas kinematics.

Zoom In Zoom Out Reset image size

Since the averages used in these profiles are mass weighted, the hot low-density interarm regions do not contribute significantly. As can be seen by inspecting Figure 2, the interarm gas has sound speeds of order $100\;\mathrm{km}\ {{\rm{s}}}^{-1}$, but as we show in Figure 9, a typical gas parcel (in a mass-weighted sense) has an effective sound speed of order $20\;\mathrm{km}\ {{\rm{s}}}^{-1}.$ We make use of mass-weighted averages so that the measured effective sound speed profiles are representative of the kinetic energy content of the gaseous disk.

Our galaxy models look quite different from real galaxies when we look at the anisotropy in the velocity dispersion components. To investigate this, in the middle panel of Figure 9 we present the time-averaged, azimuthally averaged ratio of the vertical velocity dispersion to the in-plane velocity dispersion in the fiducial galaxy model, $\sqrt{2}{\sigma }_{v,z}/{\sigma }_{v,d}.$ We include the factor of $\sqrt{2}$ in the numerator so that a value of unity would correspond to a flow field with isotropic turbulence, in which the dispersion in the two in-plane components of the velocity field account for exactly twice as much kinetic energy as the one out-of-plane component. In our model galaxies, the we find that this ratio is approximately 0.5 over the bulk of the disk, indicating that in-plane motions contribute substantially more to the turbulent kinetic energy. This can also be seen in the thickness of our disks, where we typically find that the scale height is of order one or two cell spacings on the maximum refinement level. The extreme two-dimensional (2D) disks we form are likely an artifact of the lack of feedback in our simulations.

Finally, we show in the bottom panel of Figure 9 how the effective sound speed we measure in the inner regions of the gaseous disk is primarily due to bulk turbulent motions. The small sound speed we measure is typical of the dense gas in our simulation, which cools to temperatures of $\sim$ 200–300 K due to metal line cooling. The large velocity dispersions we measure indicates that the dense gas participates in large-scale turbulent motions.

Our simulations show, in agreement with previous numerical studies by Agertz et al. (2009, 2015) and Bournaud et al. (2010), that gravitational instability is capable of driving turbulence and stabilizing a disk galaxy at $Q\gtrsim 1$ even in the absence of any additional energy input from star formation. This state is characterized by a velocity dispersion that decreases somewhat with radius from 20 km s⁻¹ near galactic centers to 10 km s⁻¹ at large radii, in good agreement with observed gas velocity dispersions in nearby disk galaxies (Tamburro et al. 2009; Ianjamasimanana et al. 2015). Indeed, the only possibly significant disagreement between the kinematic behavior found in our simulations and that observed in nature is that our galaxies' velocity dispersions are anisotropic by a factor of ∼2, leading to a scale height that is too small. Nonetheless, our results suggest that, even if stellar feedback is required to explain the vertical velocity distribution of the ISM, it need not be the dominant energy source for all turbulent motions. Gravity alone can generate the required in-plane motion.

It is interesting to consider our simulations in light of the arguments commonly made that star formation feedback is the dominant driver of turbulence in galaxies. The primary observational argument for feedback driving turbulence is that there is a correlation between star formation rate and velocity dispersion, both within galaxies (Tamburro et al. 2009) and from galaxy to galaxy (Green et al. 2010). However, our simulations would also display precisely such a correlation: the regions of our galaxy that have the highest velocity dispersion are also the most strongly star-forming, and our models with the highest gas fractions show both the highest star formation rates and the highest velocity dispersions. The point is not that feedback cannot drive turbulence, simply that a correlation between star formation and turbulence does not necessarily imply a causal relationship between the two. A correlation of this sort can be produced even when there is no feedback.

The fact that galaxies can maintain $Q\gtrsim 1$ and continuously drive strong turbulence without star formation feedback calls into question analytic models in which galaxies' star formation rates are set by the need to maintain Q = 1 (e.g., Thompson et al. 2005; Faucher-Giguère et al. 2013). We find that galaxies can maintain $Q\gtrsim 1$ regardless of the level of star formation feedback. Feedback is still needed to produce star formation rates in agreement with observation and drive out-of-plane motions, but the condition that determines the level of feedback and star formation appears to be completely decoupled from the need to maintain $Q\gtrsim 1.$

The primary result of our simulations is that gravitational instability-driven turbulence is capable of inducing significant bulk mass flows in galaxies, leading to a migration of mass inward from the passive outer regions of disks toward their actively star-forming centers. The inward mass transport rates we measure in our simulated galaxies are comparable to the star formation rates of typical L_* galaxies. While there has been a great deal of work done on such inward migration in the context of the observed giant clumps in $z\sim 2$ galaxies (e.g., Genzel et al. 2008; Bournaud & Elmegreen 2009; Cresci et al. 2009; Dekel et al. 2009b), our simulations show that a completely analogous phenomenon can operate for Milky Way-like galaxies at z = 0.

An important implication of galaxies' ability to transport gas inward at rates comparable to their star formation rates is that their available fuel for star formation is their full gas reservoir, not simply the material in the actively star-forming inner disk. This is significant because, while gas depletion times are much less than the Hubble time only considering inner disk gas, the same is not true for many galaxies at z = 0 if we consider the full gas reservoir. For a volume-limited sample of local star-forming galaxies with stellar masses in the range $\mathrm{log}({M}_{*}/{M}_{\odot })=10-11.5$, Saintonge et al. (2011) find typical H i depletion times of $\sim 3\;{\rm{Gyr}}$, with no strong dependence on stellar mass, and ${{\rm{H}}}_{2}$ depletion times that range from $\sim 0.5-3\;{\rm{Gyr}}$ from the lowest to the highest stellar masses in the sample. If galaxies have access to their full extended H i reservoirs to fuel star formation, then their total depletion times are simply the sum of these, implying that galaxies with stellar masses above $\sim {10}^{11}$ ${M}_{\odot }$ have depletion times of ∼6 Gyr. This is still less than the Hubble time, but not by much. Galaxies in this mass range could receive no new gas supply after $z\sim 1$ and still fuel all their present star formation, particularly once the contribution from stellar recycling is included. Additionally, although galaxies with masses in the range $\mathrm{log}({M}_{*}/{M}_{\odot })=10-11$ have shorter molecular gas depletion times, the atomic gas depletion times are also $\sim 3\;\mathrm{Gyr}.$ If these galaxies are also able to transport atomic gas from their outskirts into the active star-forming centers, they could supply all of their star formation without accreting any gas at all since $z\sim 0.5.$

We are therefore forced to conclude that strong equilibrium between gas inflow, star formation, and outflows may, by z = 0, exist only for dwarf galaxies with ${M}_{*}\lesssim {10}^{9.5}\;{M}_{\odot }$, since dwarfs are able to eject mass out of their star forming regions via high mass-loading factor winds. This means the mass range ${M}_{*}\sim {10}^{9.5}-{10}^{10.5}$ ${M}_{\odot }$, which corresponds to where galaxies are most efficient at converting baryons to stars (Behroozi et al. 2013), may not necessarily be in a state of near-instantaneous equilibrium between gas accretion, star formation, and outflows. Galaxies with masses above ${M}_{*}\sim {10}^{9.5}\;{M}_{\odot }$ may be far from equilibrium, contradicting one of the central assumptions of "bathtub" models of galaxy formation (e.g., Bouché et al. 2010; Lilly et al. 2013; Forbes et al. 2014b; Mitra et al. 2015).

Rather than analyzing the AMR data structures directly, we instead perform the bulk of our analysis using data interpolated onto uniform resolution meshes with 20 pc cell spacings encompassing the galactic disk out to a radius of 20 kpc and to a height of 1.25 kpc above and below the disk. This choice significantly eases the implementation of the various analysis tasks we want to perform on the data without introducing a significant amount of error. Since the disk is very thin, a vertical extent of only 1.25 kpc comfortably encloses the gas and stellar disk at all times.

To interpolate the gas data onto a uniform resolution grid, we make use of the smoothed covering grid object in yt. This operation uses cascading trilinear interpolation to represent an AMR data set at a uniform resolution. The fields representing gas density, the velocity vector components, and the gas thermal energy were extracted from the raw simulation outputs. Since the gravitational potential is not normally written to disk, we made use of the -g command-line option of the Enzo code to solve the Poisson equation using the gravity in Enzo during post-processing.

So that we can apply the same analysis tasks we use for gas fields defined on Enzo's AMR mesh to N-body star particles in our data sets, we opt to analyze stellar fields by depositing the star particle data onto the uniform resolution grids slabs. In all cases, we use cloud-in-cell interpolation onto grids with the same shape and resolution as those used for the gas data. Our choice of uniform resolution slabs leads to issues with data sparseness in the outskirts of our simulated galaxies but is a good match for the density of particle data within a radius of 10 kpc.

Because our simulations use a live stellar and dark matter halo, the rotation curve of our galaxy is not fixed, and must instead be computed self-consistently from the simulation outputs. Following Binney & Tremaine (2008) and Shu (1992), we note that for an axisymmetric system with gravitational potential field ${\rm{\Phi }}(R,z)$, we can define the circular frequency

$$\begin{eqnarray}&&{{\rm{\Omega }}}^{2}(R)=\displaystyle \frac{1}{R}{\left(\displaystyle \frac{\partial {\rm{\Phi }}}{\partial R}\right)}_{(R,0)},\end{eqnarray} \tag{ 6 }$$

and epicyclic frequency

$$\begin{eqnarray}&&{\kappa }^{2}(R)=\displaystyle \frac{2{\rm{\Omega }}}{R}\left(2R{\rm{\Omega }}+{R}^{2}\displaystyle \frac{d{\rm{\Omega }}}{{dR}}\right).\end{eqnarray} \tag{ 7 }$$

This implicitly assumes the disk is thin so we can infer the circular velocity curve by only considering the gravitational potential in the midplane.

To measure the rotation frequency, we extract the gravitational potential in a slice at the midplane of the galaxy, and evaluate the partial derivative of the potential with respect to x and y coordinates using a centered finite difference of the gravitational potential. We then form a 2D array of the gradient with respect to the cylindrical R coordinate out of the images of the x and y gradients. This results in a local estimate of the rotation frequency based on the local radial gradient in the gravitational potential. To average over local departures from axissymmetry, we create our final estimate of the rotation frequency by fitting a spline interpolator to a binned version of our local estimate of the circular frequency as a function of radius. This produces a binned 20 pc resolution estimate of the rotation frequency as a function only of cylindrical radius. Finally, we calculate the circular velocity via

$$\begin{eqnarray}&&{v}_{c}=R{\rm{\Omega }}.\end{eqnarray} \tag{ 8 }$$

Several different quantities we are concerned with are defined in terms of the projection of our simulation data. For a 3D density field $\rho (x,y,z)$ we can define the surface density,

$$\begin{eqnarray}&&{\rm{\Sigma }}(x,y)={\displaystyle \int }_{-\infty }^{\infty }\rho \;{dz}.\end{eqnarray} \tag{ 9 }$$

Here ρ represents the mass density of gas or stars, which we denote as ${\rho }_{{\rm{gas}}}$ and ${\rho }_{*}$ below. Both quantities are defined on the uniform resolution grid slabs discussed in Appendix A.1.

Discretizing ρ into a uniform resolution 3D array, which we denote as ${\rho }_{{ijk}}$, the continuous definition of the surface density reduces to

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{{ij}}=\displaystyle \sum _{k=0}^{{N}_{z}}{\rho }_{{ijk}}{\rm{\Delta }}z,\end{eqnarray} \tag{ 10 }$$

where ${\rm{\Delta }}z$ is the cell spacing the vertical direction. Since our grid slabs only include data within 1.25 kpc of the disk midplane, we are implicitly assuming that gas well off the midplane does not contribute significantly to the surface density.

We introduce the notation, ${\langle q\rangle }_{z}$, to represent the discrete mass-weighted projection operator. The result is a uniform resolution 2D array, where the $i,j$ resolution element can be found by computing,

$$\begin{eqnarray}&&{\langle q\rangle }_{z,{ij}}=\displaystyle \frac{1}{{{\rm{\Sigma }}}_{{ij}}}\displaystyle \sum _{k=0}^{{N}_{z}}{\rho }_{{ijk}}{q}_{{ijk}}{\rm{\Delta }}z.\end{eqnarray} \tag{ 11 }$$

This notation allows us to write several of the definitions below in a compact form.

Since turbulent motions are a significant component of the energy budget in the ISM, we would like to directly measure the turbulent kinetic energy in our simulations. Since our simulations are a discretized representation of a continuous underlying system, we can only estimate the velocity dispersion by comparing velocity values in cells contained within a moving window several cells across. So that variations in the rotation curve over the extent of the window do not tend to inflate the measured velocity dispersions, we define the velocity of streaming motions

$$\begin{eqnarray}&&\delta {\boldsymbol{v}}={\boldsymbol{v}}-{v}_{c}(R)\hat{\phi }.\end{eqnarray} \tag{ 12 }$$

To simplify notation in our definition of the velocity dispersion, we define the discrete operator

$$\begin{eqnarray}&&{C}_{{ijk}}(x)=\displaystyle \sum _{{k}^{\prime }=k-2}^{k+2}\displaystyle \sum _{{j}^{\prime }=j-2}^{j+2}\displaystyle \sum _{{i}^{\prime }=i-2}^{i+2}{x}_{{i}^{\prime }{j}^{\prime }{k}^{\prime }},\end{eqnarray} \tag{ 13 }$$

where x is a field defined on the uniform resolution mesh. This is equivalent to a discrete convolution of x using a cubical top-hat kernel function. The top-hat kernel includes 5 cells, so the spatial scale of the convolution is 100 pc. Using this notation, our estimate of the velocity dispersion due to motions in the x coordinate, $\delta {v}_{x}$, can be written

$$\begin{eqnarray}&&{\sigma }_{x,{ijk}}=\displaystyle \frac{\sqrt{C(\rho )C(\rho \;{\delta {v}_{x}}^{2})-C{(\rho \;\delta {v}_{x})}^{2}}}{C(\rho )},\end{eqnarray} \tag{ 14 }$$

where $\delta {v}_{x}$ is the x component of the streaming velocity vector $\delta {\bf{v}}.$ The expressions for the y and z velocity dispersions are analogous. We write the variance in terms of the convolution operator since this form lends itself to faster one-pass parallel reduction. This result can be derived from the definition of the weighted variance for a sampled quantity with measurements x_i and weights w_i, and weighted mean ${\bar{x}}_{w}$, ${\sigma }^{2}={\displaystyle \sum }_{i=1}^{N}{w}_{i}{\left({x}_{i}-{\bar{x}}_{w}\right)}^{2}/{\displaystyle \sum }_{i=1}^{N}{w}_{i}$ by expanding the squared term and substituting for the definition of the weighted mean.

We use the same definition for stellar velocity dispersions, but rather than calculating the dispersion on a cell-by-cell basis, we instead iterate over particles, performing a running dispersion calculation by depositing particle densities and velocities into accumulating arrays. This uses more memory, but is substantially faster for the particle counts in our simulations. Like the gas velocity dispersions, the stellar velocity dispersion is computed on a spatial scale of 100 pc.

Finally, we calculate the projected velocity dispersions by performing a weighted projection of the turbulent kinetic energy density:

$$\begin{eqnarray}&&{\sigma }_{v}=\sqrt{{\langle {\sigma }_{x}^{2}+{\sigma }_{y}^{2}+{\sigma }_{z}^{2}\rangle }_{z}}.\end{eqnarray} \tag{ 15 }$$

In addition, to estimate the relative contribution of in-plane and out-of-plane turbulent motions, we separately define the in-plane velocity dispersion

$$\begin{eqnarray}&&{\sigma }_{v,d}=\sqrt{{\langle {\sigma }_{x}^{2}+{\sigma }_{y}^{2}\rangle }_{z}},\end{eqnarray} \tag{ 16 }$$

and the out-of-plane velocity dispersion

$$\begin{eqnarray}&&{\sigma }_{v,z}=\sqrt{{\langle {\sigma }_{z}^{2}\rangle }_{z}}.\end{eqnarray} \tag{ 17 }$$

We can define an effective sound speed that takes into account both gas pressure and turbulent pressure,

$$\begin{eqnarray}&&{c}_{{\rm{eff}}}=\sqrt{{\sigma }_{v}^{2}+{c}_{s}^{2}}\end{eqnarray} \tag{ 18 }$$

For collisionless fluids, ${c}_{{\rm{eff}}}={\sigma }_{v}$, but for the gaseous component we must calculate the sound speed. In practice, we do this in terms of the weighted projection of the thermal energy density e,

$$\begin{eqnarray}&&{c}_{s}=\sqrt{\gamma (\gamma -1){\langle e\rangle }_{z}},\end{eqnarray} \tag{ 19 }$$

where $\gamma =5/3$ is the adiabatic index.

Using the derived data we introduced above, we can compute the Toomre Q parameter for both the gas,

$$\begin{eqnarray}&&{Q}_{{\rm{gas}}}=\displaystyle \frac{{c}_{{\rm{eff}}}\kappa }{\pi G{{\rm{\Sigma }}}_{{\rm{gas}}}}\end{eqnarray} \tag{ 20 }$$

and stars,

$$\begin{eqnarray}&&{Q}_{*}=\displaystyle \frac{{\sigma }_{v,*}\kappa }{3.36G{{\rm{\Sigma }}}_{*}}.\end{eqnarray} \tag{ 21 }$$

We also attempt to estimate the degree of gravitational instability due to both the gas and stars. We make use of the formula derived by Romeo & Wiegert (2011) that takes into account the separate contribution of the gas and stars as well as their finite thickness. In this formalism,

$$\begin{eqnarray}\displaystyle \frac{1}{{Q}_{{\rm{total}}}}=\left\{\begin{array}{lr}\displaystyle \frac{W}{{T}_{*}{Q}_{*}}+\displaystyle \frac{1}{{T}_{{\rm{gas}}}{Q}_{{\rm{gas}}}} & :{T}_{*}{Q}_{*}\geqslant {T}_{{\rm{gas}}}{Q}_{{\rm{gas}}}\\ \displaystyle \frac{1}{{T}_{*}{Q}_{*}}+\displaystyle \frac{W}{{T}_{{\rm{gas}}}{Q}_{{\rm{gas}}}} & :{T}_{*}{Q}_{*}\lt {T}_{{\rm{gas}}}{Q}_{{\rm{gas}}}\end{array}\right.,\end{eqnarray} \tag{ 22 }$$

where $T=0.8+0.7{\sigma }_{z}/{\sigma }_{d}$, and $W=2{\sigma }_{*}{\sigma }_{{\rm{gas}}}/({\sigma }_{*}^{2}+{\sigma }_{{\rm{gas}}}^{2}).$

For a fluid with local mass density ρ, we define the scale height h such that

$$\begin{eqnarray}&&{\displaystyle \int }_{0}^{h}\rho \;{dz}=\left(1-{e}^{-1}\right){\displaystyle \int }_{0}^{\infty }\rho \;{dz}.\end{eqnarray} \tag{ 23 }$$

Here the z coordinate is perpendicular to the disk and centered on the disk midplane. Since the fluid distribution is not necessarily symmetric in z, the z coordinate of the midplane is not constant, and the scale height above and below the midplane might not be equal. In practice, we take the true scale height h to be the arithmetic mean of the scale heights measured above and below the midplane.

We generate maps of of the scale height by individually processing z-aligned pencil stacks of cells. We calculate a running sum of the gas mass along each pencil stack and note the z locations where the surface density equals one of ${(2e)}^{-1}{\rm{\Sigma }}$, ${\rm{\Sigma }}/2$, and $(1.0-{(2e)}^{-1}){\rm{\Sigma }}.$ We use linear interpolation in z to estimate the intra-cell locations where the running sum exceeds each critical value. Finally, we calculate a single scale height estimate by averaging the "top" and "bottom" estimates.

For a cylindrical test volume V of radius R, the flux of mass across the surface of the cylinder is

$$\begin{eqnarray}&&\dot{M}=R{\displaystyle \int }_{0}^{2\pi }{\displaystyle \int }_{-\infty }^{\infty }\rho {v}_{r}\;{dz}\;d\theta .\end{eqnarray} \tag{ 24 }$$

In practice, we calculate this quantity using a discrete approximation based on our interpolated data defined on uniform resolution grid slabs. Here we can take advantage of the cylindrical symmetry of the problem to substantially reduce the computational cost of this calculation. If we consider a single x–y slice through our grid slab, it is a straightforward geometric problem to find the set of cells in this slice that the cylindrical test volume intersects.

If we define a one-dimensional index $l\in [0,N]$ where N is the total number of cells in the x–y plane that intersect with the circular test region, we can approximate

$$\begin{eqnarray}&&\dot{M}\approx \displaystyle \sum _{k=0}^{{k}_{{\rm{max}}}}\displaystyle \sum _{l=0}^{N}{\rho }_{{lk}}{v}_{r,{lk}}{\rm{\Delta }}{\theta }_{l}{\rm{\Delta }}z.\end{eqnarray} \tag{ 25 }$$

Here l maps to a set of unique coordinates in the 2D slice of the 3D grid slab and k indexes along the z direction in the grid slab. One can think of l as a map to a single pair of $i,j$ indices for the x and y rows in our grid slice. The problem of calculating the radial mass flux reduces to calculating the angles subtended by each cell ${\rm{\Delta }}{\theta }_{l}.$ Since we ignore gas well off the midplane, this implicitly assumes that the mass transport at large heights above and below the midplane is small.

Doing this at many different radii allows up to map how the radial mass flux varies with galactocentric radius. In practice, we choose a set of 1000 nested test cylinders, binned evenly in radius with a spacing of 40 pc.

We are interested in the time-average behavior of azimuthally averaged quantities, both to compare to 1D models of galaxy formation and also to capture the large-scale behavior of our simulations over long periods of time. As we show in the main text, the internal structure of our simulated galaxies is initially dominated by transient disturbances. After several galactic rotation periods, the disks settle down and reach a quasi-equilibrium state in which the global structure is more or less static. Since the galactic rotation period is an increasing function of radius, the places where we expect the disk to be settled corresponds to a wedge-shaped region of radius-time phase space.

For any azimuthally averaged quantity $x(r,t)$, we define its time average, a function of radius only, as

$$\begin{eqnarray}&&{\langle x(r,t)\rangle }_{t}=\displaystyle \frac{1}{{t}_{{\rm{max}}}-({t}_{0}+{ar})}{\displaystyle \int }_{{t}_{0}+{ar}}^{{t}_{{\rm{max}}}}x(r,t){dt},\end{eqnarray} \tag{ 26 }$$

where ${t}_{{\rm{max}}}$ is the simulation time at the end of the simulation (600 Myr for all three simulations presented here) and a and t₀ are constants that define the shape of the "averaging wedge" as denoted by the blue lines in Figure 6. In practice we determine these constants "by eye" based on inspection of the simulation results.